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METHOD AND SYSTEM FOR TRACKING MULTIPLE REGIONAL 
OBJECTS BY MULTI-DIMENSIONAL RELAXATION 



[FIELD OF THE INVENTION 



The invention relates generally to computerized techniques 
for processing data obtained from radiation reflections lASed to 
track multiple discrete object. ^v*£> 



This application is a continuation of U.S. Patent 
Application Serial No. 09/312,036 filed May 14, 1999 which is 
a continuation of U.S. Patent Application Serial No. 08/682,904 
filed July 16, 1996, now U.S. Patent 5,959,574 which is a 
continuation-in-part of U.S. Patent Application Serial No. 
08/404,024, filed March 14, 1995 and issued July 16, 1996 as 
U.S. Patent No. 5,537,119, which is a continuation-in-part of 
U.S. Patent Application Serial No. 08/171,327 filed December 
21, 1993, now U.S. Patent No. 5,406,289. 



a . Field of the Invention 

The invention relates generally to computerized techniques 
for processing data obtained from radar to track multiple 
discrete objects. 



1 CROSS REFERENCE TO RELATED APPLICATIONS 




BACKGROUND OF THE INVENTION 
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There are many situations where the courses of multiple 
objects in a region must be tracked. Typically, radar is used 
to scan the region and generate discrete images or "snapshots" 
based on sets of returns or observations. In some types of 
5 tracking systems, all the returns from any one object are 
represented in an image as a single point unrelated to the 
shape or size of the objects. "Tracking" is the process of 
identifying a sequence of points from a respective sequence of 
the images that represents the motion of an object. The 

10 tracking problem is difficult when there are multiple closely 
spaced objects because the objects can change [ their] speed and 
direction rapidly and move into and out of the line of sight 
for other objects. The problem is exacerbated because each set 
of returns may result from noise as well as echoes from the 

15 actual objects. The returns resulting from the noise are also 
called false positives. Likewise, the radar will not detect 
all echoes from the actual objects and this phenomena is called 
a false negative or "missed detect" error. For tracking 
airborne objects, a large distance between the radar and the 

20 objects diminishes the signal to noise ratio so the number of 
false positives and false negatives can be high. For robotic 
applications, the power of the radar is low and as a result, 
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the signal to noise ratio can also be low and the number of 
false positives and false negatives high. 

In view of the proximity of the objects to one another, 
varied motion of the objects and false positives and false 
5 negatives, multiple sequential images should be analyzed 
collectively to obtain enough information to properly assign 
the points to the proper tracks. Naturally, the larger the 
number of images that are analyzed, the greater the amount of 
information that must be processed. 

10 While identifying the track of an object, a kinematic 

model may be generated describing the [ object's] location, 
velocity and acceleration [may be generated! of the object . 
Such a model provides the means by which the [object's ] future 
motion of the object can be predicted. Based upon such a 

15 prediction, appropriate action may be initiated. For example, 
in a military application there is a need to track multiple 
enemy aircraft or missiles in a region to predict their 
objective, plan responses and intercept them. Alternatively, 
in a commercial air traffic control application there is a need 

20 to track multiple commercial aircraft around an .airport to 
predict their future courses and avoid collision. Further, [ 
in] these and other applications, such as robotic 
applications, may use radar, sonar, infrared or other object 
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detecting radiation bandwidths for tracking objects. In 
particular, in robotic applications reflected radiation can be 
used to track a single object which moves relative to the robot 
(or vice versa) so the robot can work on the object. 
5 Consider the very simple example of two objects being 

tracked and no false positives or false negatives. The radar, 
after scanning at time t lf reports objects at two locations in 
a first observation set. That is, it returns a set of two 
observations {o llf o 12 ) . At time t 2 it returns a similar set of 

10 two observations {o 21 , o 22 ) from a second observation set. 
Suppose from prior processing that track data for two tracks T 1 
and T 2 includes the locations at t 0 of two objects. Track T x 
may be extended through the points in the two sets of 
observations in any of four ways, as may track T 2 . The possible 

15 extensions of T 2 can be described as: {T lf o n , o 21 } , {T lf o llf 
o 22 ) , {T lf o 12 , o 21 } and { T 1 , o 12r o 22 ) . Tracks can likewise be 
extended from T 2 in four possible waySj. including [,] { T 2 , o 12 , 
o 21 } . Figure 1 illustrates these five (out of eight) possible 
tracks (to simplify the problem for purposes of explanation) . 

20 The five track extensions are labeled h n , h 12 , h 13 , h 14r and h 21 
wherein h n is derived from {T lf o n , o 21 ) f h 12 is derived from 
{T lf o llf o 22 ) , h 13 is derived from {T lr o 12 , o 21 ] , h 14 is derived 
from {T lf o 12 , o 22 ) f and h 21 is derived from { T 2 , o n , o 21 ) . The 
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problem of determining which such track extensions are the most 
likely or optimal is hereinafter known as the assignment 
problem. 

It is known from prior art to determine a figure of merit 
5 or cost for assigning each of the points in the images to a 
track. The figure of merit or cost is based on the likelihood 
that the point is actually part of the track. For example, the 
figure of merit or cost may be based on the distance from the 
point to an extrapolation of the track. Figure 1 illustrates 

10 costs 5 21 and 5 22 for hypothetical extension h 21 and modeled 
target characteristics. The function to calculate the cost 
will normally incorporate detailed characteristics of the 
sensor, such as probability of measurement error, and track 
characteristics, such as likelihood of track maneuver. 

15 Figure 2 illustrates a two by two by two matrix, c, that 

contains the costs for each point in relation to each possible 
track. The cost matrix is indexed along one axis by the track 
number, along another axis by the image number and along the 
third axis by a point number. Thus, each position in the cost 

20 matrix lists the cost for a unique combination of points and a 
track, one point from each image. Figure 2 also illustrates a 
{0, 1} assignment matrix, z, which is defined with the same 
dimensions as the cost matrix. Setting a position in the 
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assignment matrix to "one" means that the equivalent position 
in the cost matrix is selected into the solution. The 
illustrated solution matrix selects the {h l4 , h 21 ) solution 
previously described. Note that for the above example of two 
5 tracks and two snapshots, the resulting cost and assignment 
matrices are three [ ] -dimensional . As used in this patent 
application, the term "dimension" means the number of axes in 
the cost or assignment matrix while size refers to the number 
of elements along a typical axis. The costs and assignments 

10 have been grouped in matrices to facilitate computation. 

A solution to the assignment problem satisfies two 
constraints - first, the sum of the associated costs for 
assigning points to a track extension is minimized and, second, 
if no false positives or false negatives exist, then each point 

15 is assigned to one and only one track. 

When false positives exist, however, additional 
hypothetical track extensions incorporating the false positives 
will be generated. Further note that the random locations of 
false positives will, in general, not fit well with true data 

20 and such additional hypothetical track extensions will result 
in higher costs. Also note that when false negative errors 
exist, then the size of the cost matrix must grow to include 
hypothetical track extensions formulated with "gaps" (i.e., 
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data omissions where there should be legitimate observation 
data) for the false negatives. Thus, the second criteria must 
be weakened to reflect false positives not being as signed and 
also to permit the gap filler to be multiply assigned. With 
hypothetical cost calculated in this manner then the foregoing 
criteria for minimization will tend to materialize the false 
negatives and avoid the false positives. 

For a r 3-dimensionall three-dimensional problem, as is 
illustrated in Figure 1, but with N 2 (initial) tracks, N 2 
observations in scan 1, N 3 observations in scan 2, false 
positives and negatives assumed, the assignment problem can be 
formulated as: 

N, N~ «L 

(a) Minimize: NL E E C W , Z W, 

i^O i 2 =0 i 3 =0 123 123 

(b) Subject to: [] X, Lj z t ± ± =1 ' i i = 1 ' ■ • • N i-^- 



i 2 =i 1,-1 ^ 



(C) UX, 22 Z iii^ 1 ' i 2 = 

i 3 =l 



2-c 



(d) []22 12 tin* 1 * i 3 = l,...N^_ 

i 1= l i 2 =l 1 2 3 

(e) n* W 3 e{0 '' 1} nv z w/ 



[ [1.0] ] 
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where "c" is the cost and "z" is a point or observation 
assignment, as in Figure 2. 

The minimization equation or equivalently objective 
function \ \] (0 . 1 \ . 01 1 ) (a) specifies the sum of the element by 
5 element product of the c and z matrices. The summation 
includes hypothesis representations z i 1 i 2 i 3 with observation 
number zero being the gap filler observation. Equation 
r n (0. 1 f .011 ) (b) requires that each of the tracks T lf ...,T Nl be 
extended by one and only one hypothesis. Equation 

10 r n ( 0 . 1 f . 011 ) (c) relates to each point or observation in the 
first observation set and requires that each such observation, 
except the gap filler, can only associate with one track but^ 
because of the "less than" condition^ it might not associate 
with any track. Equation HI (0. 1 \ . 011 ) (d) is like 

15 r fl ( 0 . 1 f . 01 1 ) (c) except that it is applicable to the second 
observation set. Equation f [1 (0 . 11". 011) (e) requires that 
elements of the solution matrix z be limited to the zero and 
one values. 

The only known method to solve Problem Formulation 
20 rn (0. U.011) exactly is a method called "Branch and Bound." 
This method provides a systematic ordering of the potential 
solutions so that solutions with a same partial solution are 
accessible via a branch of a tree describing all possible 
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solutions whereby the cost of unexamined solutions on a branch 
are incrementally developed as the cost for other solutions on 
the branch are determined. When the developing cost grows to 
exceed the previously known minimal cost (i.e., the bound) then 
5 enumeration of the tree branch terminates. Evaluation 
continues with a new branch. If evaluation of the cost of a 
particular branch completes, then that branch has lower cost 
than the previous bound so the new cost replaces the old bound. 
When all possible branches are evaluated or eliminated then the 

10 branch that had resulted in the last used bound is the 
solution. If we assume that N x = N 2 = N 3 = n and that branches 
typically evaluate to half there full length, then workload 
associated with "branch and bound" is proportional to 
(n ! | [ ] ! ) 2 [ ] . This workload is unsuited to real time 

15 evaluation. 

The Branch and Bound algorithm has been used in past 
research on the Traveling Salesman Problem. Messrs. Held and 
Karp showed that if the set of constraints was relaxed by a 
method of Lagrangian Multipliers (described in more detail 
20 below) then tight lower bounds could be developed in advance of 
enumerating any branch of the potential solution. By 
initiating the branch and bound algorithm with such a tight 
lower bound, significant performance improvements result in 
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that branches will typically evaluate to less than half their 
full length, 

Messrs. Frieze and Yadagar in dealing with a problem 
related to scheduling combinations of resources, as in job, 
5 worker and work site, showed that Problem Formulation 
rn ( 0 . 1 I" . 01 1 ) applied. They further described a solution 
method based upon an extension of the Lagrangian Relaxation 
previously mentioned. The two critical extensions provided by 
Messrs. Frieze and Yadagar were: (1) an iterative procedure 

10 that permitted the lower bound on the solution to be improved 
(by "hill climbing" described below) and (2) the recognition 
that when the lower bound of the relaxed problem was maximized, 
then there existed a method to recover the solution of the non- 
relaxed problem in most cases using parameters determined at 

15 the maximum. The procedures attributed to Messrs. Frieze and 
Yadagar are only applicable to the f 3-dimensionall three- 
dimensional problem posed by Problem Formulation rn ( 0 . If. 011) 
and where the cost matrix is fully populated. However, 
tracking multiple airborne objects usually requires solution of 

20 a much higher dimensional problem. 

Figures 1 and 2 illustrate an example where "look ahead" 
data from the second image improved the assignment accuracy for 
the first image. Without the look ahead, and based only upon 
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a simple nearest neighbor approach, the assignments in the 
first set would have been reversed. Problem Formulation 
I'M (O. ir.011) and the prior art only permit looking ahead one 
image. In the prior art it was known that the accuracy of 
5 assignments will improve if the process looks further ahead, 
however no practical method to optimally incorporate look ahead 
data existed. Many real radar tracking problems involve 
hundreds of tracks, thousands of observations per observation 
set and matrices with dimensions in the range of 3 to 25 

10 including many images of look ahead. 

It was also known that the data assignment problem may be 
simplified (without reducing the dimension of the assignment 
problem) by eliminating from consideration for each track those 
points which, after considering estimated limits of speed and 

15 turning ability of the objects, could not physically be part of 
the track. One such technique, denoted hereinafter the "cone 
method, " defines a cone as a continuation of each previously 
determined track with the apex of the cone at the end of the 
previously defined track. The length of the cone is based on 

20 the estimated maximum speed of the object and the size of the 
arc of the cone is based on the estimated maximum turning 
ability of the object. Thus, the cone defines a region outside 
of which no point could physically be part of the respective 
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track. For any such points outside of the cones, an infinite 
number could be put in the cost matrix and a zero could be 
preassigned in the assignment matrix. It was known for the 
tracking problem that these elements will be very common in the 
5 cost and selection matrices (so these matrices are "sparse"). 

It was also known in the prior art that one or more tracks 
which are substantially separated geographically from the other 
tracks can be separated also in the assignment problem. This 
is done by examining the distances from each point to the 

10 various possible tracks. If the distances from one set of 
points are reasonably short only in relation to one track, then 
they are assigned to that track and not further considered with 
the remainder of the points. Similarly, if a larger group of 
points can only be assigned to a few tracks, then the group is 

15 considered a different assignment problem. Because the 
complexity of assignment problems increases dramatically with 
the number of possible tracks and the total number of points in 
each matrix, this partitioning of the group of points into a 
separate assignment problem and removal of these points from 

20 the matrices for the remaining points, substantially reduces 
the complexity of the overall assignment problem. 

A previously known Multiple Hypothesis Testing (MHT) 
algorithm (see Chapter 10 of S.S. Blackman[,1. Multiple-Target 
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Tracking with Radar Applications [ , Chapter !(),]_,_ Artech House^ 
Norwoods MA, 1986^_) related to formulation and scoring. The 
MHT procedure describes how to formulate the sparse set of all 
reasonable extension hypothesis (for Figure 1 the set 
5 {h u . . .h 24 }) and how to calculate a cost of the hypothesis {T ± , 
°2k) based upon the previously calculated cost for 
hypothesis {T if o 2j *} . The experience with the MHT algorithm, 
known in the prior art, is the basis for the assertion that 
look ahead through k sets of observations results in improved 

10 assignment of observations from the first set to the track. 

In theory, the MHT procedure uses the extendable costing 
procedure to defer assignment decision until the accumulated 
evidence supporting the assignment becomes overwhelming. When 
it makes the assignment decision^, it then eliminates all 

15 potential assignments invalidated by the decision in a process 
called "pruning the tree." In practice, the MHT hypothesis 
tree is limited to a fixed number of generations and the 
overwhelming evidence rule is replaced by a most likely and 
feasible rule. This rule considers each track independently of 

20 others and is therefore a local decision rule. 

A general object of the present invention is to provide an 
efficient and accurate process for assigning each point object 
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in a region from multiple images to a proper track and then 
taking an action based upon the assignments. 

A more specific object of the present invention is to 
provide a technique of the foregoing type which determines the 
5 solution of a ^--dimensional assignment problem where "Jc" is 
greater than or equal to three. 

SUMMARY OF THE INVENTION 
The present invention relates to a method and apparatus 

10 for tracking objects. In particular, the present invention 
tracks movement or trajectories of obj.ects by analyzing 
radiation reflected from the objects, the invention being 
especially useful for real-time tracking in noisy environments. 
In providing such a tracking capability, a region 

15 containing the objects is repeatedly scanned to generate a 
multiplicity of sequential images or data observation sets of 
the region. One or more points (or equivalently observations) , 
in each of the images or observation sets are detected wherein 
each such observation either corresponds to an actual location 

20 of an object or is an erroneous data point due to noise. 
Subsequently, for each observation detected, figures of merit 
or costs are determined for assigning the observation to each 
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of a plurality of previously determined tracks. *Afterwards, 
a first optimization problem is specified which includes: 

(a) a first objective function for relating the above 
mentioned costs to potential track extensions through the 

5 detected observations (or simply observations) ; and 

(b) a first collection of constraint sets wherein each 
constraint set includes constraints to be satisfied by the 
observations in a particular scan to which the constraint 
set is related. In general, there is a constraint for 

10 each observation of the scan, wherein the constraint 

indicates the number of track extensions to which the 
observation may belong. 

In particular, the first optimization problem is 
formulated, generated or defined as an Af-dimensional assignment 

15 problem wherein there are M constraint sets in the first 
collection of constraint sets (i.e., there are M scans being 
examined) and the first objective function minimizes a total 
cost for assigning observations to various track extensions 
wherein terms are included in the cost, such that the terms 

20 have the figures of merit or costs for hypothesized 
combinations of assignments of the observations to the tracks. 
Subsequently, the formulated M-dimensional assignment problem 
is solved by reducing the complexity of the problem by 
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generating one or more optimization problems each having a 
lower dimension and then solving each lower dimension 
optimization problem. That is, the M-dimensional assignment 
problem is solved by solving a plurality of optimization 
5 problems each having a lower number of constraint sets. 

The reduction of the M-dimensional assignment problem to 
a lower dimensioned problem is accomplished by relaxing the 
constraints on the points of one or more scans thereby 
permitting these points to be assigned to more than one track 

10 extension. In relaxing the constraints, terms having penalty 
factors are added into the objective function thereby 
increasing the total cost of an assignment when one or more 
points are assigned to more than one track. Thus, the 
reduction in complexity by this relaxation process is 

15 iteratively repeated until a sufficiently low dimension is 
attained such that the lower dimensional problem may be solved 
directly by known techniques. 

In one embodiment of the invention, each £:-dimensional 
assignment problem f2<k^M, 1 (2 < k < M) is iteratively reduced 

20 to a fk-11 (k-1) - dimensional problem until a [2- 
dimensionall two-dimensional problem is specified or 

formulated. Subsequently, the r2-dimensionall two-dimensional 
problem formulated is solved directly and a "recovery" 
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technique is used to iteratively recover an optimal or near- 
optimal solution to each ic-dimensional problem from a derived 
(k-1) [ dimensional! -dimensional problem r k=2 1 k= 2 ,3,4, . , >M. 

In performing each recovery step (of obtaining a solution 
5 to a ^-dimensional problem using a solution to a (k-1)- 
dimensional problem) ^_ an auxiliary function [,] is utilized. In 
particular, to recover an optimal or near-optimal solution to 

a A:-dimensional problem, an auxiliary _function, _ff J j t _ 

Ck=A, 5, . . . , AQ_, is specified and a region or domain is 

10 determined wherein this function is maximized, whereby values 
of the region determine the penalty factors of the (k-1) - 
dimensional problem such that another r2-dimensional1 two- 
dimensional problem can be formulated which determines a 
solution to the £:-dimensional problem using the penalty factors 

15 of the (k-1) -dimensional problem. 

Each Auxiliary function ¥ kt k=3, . . . ,M-1, is a function of 
both lower dimensional problems , penalty factors^ and a 
solution at the dimension k at which the penalized cost 
function is solved directly (typically a r2-dimensional1 two- 

20 dimensional problem) . Further, in determining, for auxiliary 
function W kf a peak region, a gradient of the auxiliary 
function is determined, and utilized to identify the peak 
region. Thus, gradients are used for each of the approximation 
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functions W 3f W 4f . . . , W M _j which are sequentially formulated and 
used in determining penalty factors [penalty factors ] until 
\ tm- 1 1 M-l is used in determining the penalty factors for the 
[M-l dimensional! (M-l ) -dimensional problem. Subsequently, 
5 once the M-dimensional problem is solved (using a [2- 
dimensionall two-dimensional problem to go from an (M-l) [ 
dimensional] -dimensional solution to an M-dimensional 
solution) , one or more of the following actions are taken 
based on the track assignments: sending a warning to aircraft 

10 or a ground or sea facility, controlling air traffic, 
controlling anti-aircraft or anti-missile equipment, taking 
evasive action, working on one of the objects. 

According to one feature of this first embodiment of the 
present invention, the following steps are also performed, 

15 before the step of defining the auxiliary function. A 
preliminary auxiliary function is defined for each of the lower 
dimensional problems having a dimension equal or one greater 
than the dimension at which the penalized cost function is 
solved directly. The preliminary auxiliary function is a 

20 function of lower order penalty values and a solution at the 
dimension at which the penalized cost function was solved 
directly. In determining a gradient of the preliminary 
auxiliary function, step in the direction of the gradient to 
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identify a peak region of the preliminary auxiliary function 
and determine penalty factors at the peak region. Iteratively 
repeat the defining, gradient determining, stepping and peak 
determining steps to define auxiliary functions at successively 
5 higher dimensions until the auxiliary function [at 6-18 (k-1) 
1 dimension lk-1) is determined. In an alternative second 
embodiment of the present invention, instead of reducing the 
dimentionality of the M-dimensional assignment problem by a 
single dimension at a time, a plurality of dimensions are 

10 relaxed simultaneously. This new strategy has the advantage 
that when the M-dimensional problem is relaxed directly to a [ 
2-dimensionall two-dimensional assignment problem, then all 
computations may be performed precisely without utilizing an 
auxiliary function such as W k as in the first embodiment. More 

15 particularly, the second embodiment solves the first 
optimization problem (i.e., the [M- dimensional"! M-dimensional 
assignment problem) by specifying (i.e., creating, generating, 
formulating and/or defining) a second optimization ^.problem. [ 
] The second optimization problem includes a second objective 

20 function and a second collection of constraint sets wherein: 
a) the second objective function is a combination of the 
first objective function and penalty factors or terms 
determined for incorporating ±M-m [ constraint 1 ) -constraint 
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sets of the first optimization problem into the second 
obj ective function; 
b) the constraint sets of the second collection include only 
m of the constraint sets of the first collection of 
5 constraints, 

wherein r2<m<M-11 2 < m < M -1 . Note that, once the second 
optimization problem has been specified or formulated, an 
optimal or near-optimal solution is determined and that 
solution is used in specifying (i.e, creating, generating, 

10 formulating and/or defining) a third optimization problem of 
±M-m}_ dimensions (or equivalently constraint sets) . The third 
optimization problem is subsequently solved by decomposing it 
using the same procedure of this second embodiment as was used 
to decompose the first optimization problem above. Thus, a 

15 plurality of instantiations of the third optimization problem 
are specified, each successive instantiation having a lower 
number of dimensions, until an instance of the third 
optimization problem is a two [ ] -dimensional assignment problem 
which can be solved directly. Subsequently, whenever an 

20 instance of the third optimization problem is solved, the 
solution is used to recover a solution to the instance of the 
first optimization problem from which this instance of the 
third optimization was derived. Thus, an optimal or near- 
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optimal solution to the original first optimization problem may 
be recovered through iteration of the above steps. 

As mentioned above, the second embodiment of the present 
invention is especially advantageous when \m=2] m =2 , since in 
5 this case all computations may be performed precisely and 
without utilizing auxiliary functions. 

Other features and benefits of the present invention will 
become apparent from the detailed description with the 
accompanying drawings contained hereinafter. 

10 

BRIEF DESCRIPTION OF THE FIGURES 
Figure 1 is a graph of images or data sets generated by a 
scan of a region and possible tracks within the images or data 
sets according to the prior art. 
15 Figure 2 illustrates cost and assignment matrices for the 

data sets of Figure 1 according to the prior art. 

Figure 3 is a block diagram of the present invention. 
Figure 4 is a flow chart of a process according to the 
prior art for solving a r 3 -dimensional 1 three-dimensional 
20 assignment problem. 

Figure 5 [I]is a flow chart of a process according to the 
present invention for solving a Jc-dimensional assignment 
problem where "Jr" is greater than or equal to 3. 
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Figure 6 is a graph of various functions used to explain 
the present invention. 

Figure 7 is another block diagram of the present invention 
for solving the ^-dimensional assignment problem where "k" is 
5 greater than or equal to three ( 3) . 

Figure 8 is a flowchart describing the procedure for 
solving an n-dimensional assignment problem according to the 
second embodiment of the invention. 

10 DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 

Referring now to the other figures in detail wherein like 
reference numerals indicate like elements throughout the 
several views, Figure 3 illustrates a system generally 
designated 100 for implementing the present invention. System 

15 100 comprises, for example, a radar station 102 (note sonar, 
microwave, infrared and other radiation bandwidths are also 
contemplated) for scanning a region which may be, for example, 
an aerial region (in aerial surveillance applications) or a 
work region (in robotic applications) and generating signals 

20 indicating locations of objects within the region. The signals 
are input to a converter 104 which converts the signals to data 
points or observations in which each object (or false positive) 
is represented by a single point. The output of the converter 
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is input to and readable by a computer 10 6. As described in 
more detail below, the computer 106 assigns the points to 
respective tracks, and then displays the tracks and 
extrapolations of the tracks on a monitor 110. Also, the 
5 computer 106 determines an appropriate action to take based on 
the tracks and track extensions. For example, in a commercial 
application at an airport, the computer can determine if two 
aircraft being tracked are on a collision course and if so, 
signal a transmitter 112 to warn each aircraft, or if a 

10 scheduled take-off will pose the risk of collision, delay the 
take-off. For a military application on a ship or base, the 
computer can determine subsequent coordinates of enemy aircraft 
and send the coordinates to an antiaircraft gun or missile 120 
via a communication channel 122. In a robotic application, the 

15 computer controls the robot to work on the proper object or 
portion of the object. 

The invention generates A: -dimensional matrices where k is 
the number of images or sets of observation data in the look 
ahead window plus one. Then, the invention formulates a k- 

20 dimensional assignment problem as in r I" 1 (0 . 1 f . 01 1 ) above. 

The ^-dimensional assignment problem is subsequently 
relaxed to a (Jc-1) -dimensional problem by incorporating one set 
of constraints into the objective function using a Lagrangian 
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relaxation of this set. Given a solution of the {k-1) - 
dimensional problem, a feasible solution of the ^-dimensional 
problem is then reconstructed. The (k-1 ) -dimensional problem 
is solved in a similar manner, and the process is repeated 
5 until it reaches the two-dimensional problem that can be solved 
exactly. The ideas behind the Lagrangian relaxation scheme are 
outlined next. 

Consider the integer programming problem 
Minimize v(z) = c T z 

10 

Subject to: Az < b^ (0.2) 

Bz < cL 

z i is an integer for i G Ij_ 

[[1.0.1]] 

15 where the partitioning of the constraints is natural in some 
sense. Given a multiplier vector [u>01 u > 0 , the Lagrangian 
relaxation of [ [1 . ] ±0 . [1] ] 2J_ relative to the constraints 
fBz<d] Bz < d is defined to be 

0(u) = Minimize { c T z+u T (Bz-d) } (0.3) 

20 

Subject to: [Az<bl Az < b, 

z i is an integer for i E 

[[1.0.2]] 

If the constraint set fBz^dl Bz < d is replaced by fBz=d] Bz 

25 = d, the non-negativity constraint on u is removed. 

^=c T z-f-u T (Bz-d) is the Lagrangian relative to the constraints 
[Bz^dl Bz d , and hence the name Lagrangian relaxation. The 
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following fact gives the relationship between the objective 
functions of the original and relaxed problems. 

FACT A.l. If F is an optimal solution to f fl. 1 (0. fll 12) , 
then <P(u) [ <1 < v(z~) for all fu>01 u > 0 . If an optimal solution z 
5 of [ [1. ]J_0. [2] ]3l is feasible for [ [1 . ] 10 . [2] ] 3L, then z is an 
optimal solution for [ [ 1 . ] ±0 . [ 1 ] ] 2± and <P(u) [=1= v(z) . 
This leads to the following algorithm: 

Algorithm. Construct a sequence of multipliers {u k }% =0 
converging to the solution u of Maximize {<P{u) : ru>01 u > 0 ) 
10 and a corresponding sequence of feasible solutions {z k }^ 0 of 
[ [1.110. [1] ]21 as follows: 

1. Generate [in] an initial approximation u 0 . 

2. Given u k , choose a search direction s k and a search 
distance a k so that <P{ [u k +a k s k ] u k _t& k s k ) [>] >__<&{ u k ) . Update the 

15 multiplier estimate u k by ru k+1 =u k +a k s k 1 u y+1 = u k + a k s^. 

3. Given u k+1 and a feasible solution Jik+i( u k+i) of 
[ [1 . ] J_0 . [2] ] 3J_, recover a feasible solution z k+1 {u k+1 ) of the 
integer programming problem [ [1 . ]±0. [1] ] 2J_. 

4. Check the termination criteria. If the algorithm is 
20 not finished, set k=k+l and return to Step 2. Otherwise, 

terminate the algorithm. 

If F is an optimal solution of [ [1 . ]±0 . [1] ] 2J_, then 
<P(u k ) <> <P(u) <> v(z~) <> v(z k )^ 
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Since the optimal solution z = f(u) of [ [1 • ] 10 . [2 ] ] 3JL is 

usually not a feasible solution of [ [ 1 . ] 10 . [ 1 ] ] 2_]_ for any 
choice of the multipliers, <P(u) is usually strictly less than 
v(F). The difference v(F)-<P(u) is called the "duality gap, " 
5 which can be estimated by comparing the best relaxed and 
recovered solutions computed in the course of maximizing <t>(u). 

Thus, to eliminate the "<" in the constraint formulations 
(e.g., r fl (0.1 r .01 1 ) (a) - [ []_UL_1[.0] ]l(d) ) that resulted from 
false positives, the constraint sets have been modified, as 

10 indicated above, to incorporate a gap filler in each 
observation set to account for false positives. Thus, a zero 
for an index [, ]_i p /", l<p<kl (1 < p < k) in z k ilm _ ±k in problem 
[ [1J1Q . [1] ] 4J_ below indicates that the track extension 
represented by z k ± ± includes a false positive in the p th 

15 observation set. Note that this implies that a hypothesis be 
formed incorporating an observation with k-1 gap fillers, e.g., 
z o...oi p o...or i P *0. Thus, the resulting generalization of Problem 
Formulation [ f 1 (0. IT. 011) without the "less than" complication 
within the constraints is the following Jc-Dimensional 

20 Assignment Problems in which \ k> 3 1 k z3 : 
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[ Minimize :1 Minimize []/>••• Z-/ c v z < 
[ ] Subject to: []£••■ 2 z. . = 1, [ ] i 



i 2 =o i k =o -I"-* 



(0.4) 



5^ ••• £ £ jtl Z i,...i u ^ / 

i 1= 0 1^ = 0 i J+1 =0 vo 1 

5 [ 

] [ ]for = 1,. . .Nj 

[ ] and j=2,...k 



- 1 



[ ] []£••• £ Z i i=lr [ ]i k = l,...N kJ _ 

10 [ ] []z V--* 6 {0 ' 1} ' 

[ [1.1]] 

where c and z are similarly dimensioned matrices representing 
costs and hypothetical assignments. Note, in general, for 

15 tracking problems these matrices are sparse. 

After formulating the A:-dimensional assignment problem as 
in [ [1] J_0. [1] ] 4J_, the present invention solves the resulting 
problem so as to generate the outputs required by devices 110, 
112, 122 and 130. For each observation set received from 

20 converter 104 at time t i where i=l,...,A/ 2 , the computer 
processes Oi in a batch together with the other observation sets 
°i-k+if • • • / °i and the track T ± . kf i.e., Tj is the set of all tracks 
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that have been defined up to but not including Oj . (Note, bold 
type designations refer to the vector of elements for the 
indicated time, i.e., the set of all observations in the scan 
or tracks existing at the time, etc.) The result of this 
5 processing is the new set of tracks Ti„ k+1 and a set of cost 
weighted possible solutions indicating how the tracks might 
extend to the current time t ± . At time t ±+1 the batch process 
is repeated using the newest observation set and deleting the 
oldest. Thus, there is a moving window of observation sets 

10 which is shifted forward to always incorporate the most recent 
observation set. The effect is that input observation sets are 
reused for rk-11 (k-1) cycles and then on the observation set's 
[k-th]ir^ reuse each observation within the observation set is 
integrated into a track. 

15 Figure 7 illustrates various processes implemented upon 

receipt of each observation set. Except for the addition of 
the A:-dimensional assignment solving process 300 and the 
modification to scoring process 154 to build data structures 
suitable for process 300, all processes in Figure 7 are based 

20 upon prior art. The following processes 150 and 152 extend 
previously defined tracks h^^.] based on new observations. Gate 
formulation and output process 156 determines, for each of the 
previously defined tracks, a zone wherein the track may 
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potentially extend based on limits of velocity, maneuverability 
and radar precision. One such technique to accomplish this is 
the cone method described previously. The definition of the 
zone is passed to gating process 150. When a new observation 
5 set Oi is received, the gating process 150 will match each 
member observation with the zone for each member of the 
hypothetical set h ± _ ± . After all input observations from 0 ± are 
processed^, the new hypothesis set h ± is generated by extending 
each track of the prior set of hypothetical tracks h^ either 

10 with missed detect gap fillers or with all new observation 
elements satisfying the track's zone. This is a many [ ]-to[ ]- 
many matching in that each hypothesis member can be extended to 
many new observations and each new observation can be used to 
extend many hypotheses. It, however, is not a full matching in 

15 that any hypothesis will neither be matched to all observations 
nor vice versa. It is this matching characteristic that leads 
to the sparse matrices involved in the tracking process. 
Subsequently, gating 150 forwards the new hypothesis set h ± to 
filtering process 152. Filtering process 152 determines a 

20 smooth curve for each member of hi. Such a smooth curve is more 
likely than a sharp turn from each point straight to the next 
point. Further, the filtering process 152 removes small errors 
that may occur in generating observations. Note that in 
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performing these tasks, the filtering process 152 preferably 
utilizes a minimization of a least squares test of the points 
in a track hypothesis or a Kalman [F] filtering approach. 

As noted above, the foregoing track extension process 
5 requires knowledge of a previous track. For the initial 
observations, the following gating process 158 and filtering 
process 160 determine the "previous track" based on initial 
observations. In determining the initial tracks, the points 
from the first observation set form the beginning points of all 

10 possible tracks. After observation data from the next 
observation set is received, sets of simple two [ ] -point 
straight line tracks are defined. Then, promotion, gate 
formulation, and output step 162 determines a zone in which 
future extensions are possible. Note that filtering step 160 

15 uses curve fitting techniques to smooth the track extensions 
depending upon the number of prior observations that have 
accumulated in each hypothesis. Further note that promotion, 
gate formulation and output process 162 also determines when 
sufficient observations have accumulated to form a basis for 

20 promoting the track to processes 150 and 152 as described 
above . 

The output of each of the filtering processes 152 and 160 
is a set of hypothetical track extensions. Each such extension 
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contains the hypothetical initial conditions (from the previous 
track) , the list of observations incorporated in the extension, 
and distance between each observation and the filtered track 
curve. Scoring process 154 determines the figure of merit or 
5 cost of an observation being an extension of a track. In one 
embodiment, the cost is based on the above-mentioned distance 
although the particular formula for determining the cost is not 
critical to the present invention. A preferred formula for 
determining the cost utilizes a negative log likelihood 

10 function in which the cost is the negative of the sum of: (a) 
the logs of the distances normalized by sensor standard 
deviation parameters, and (b) the log likelihoods for events 
related to [ : ] track initiation, track termination, track 
maneuver, false negatives and false positives. Note that track 

15 maneuvers are detected by comparing the previous track curve 
with the current extension. Further note that some of the 
other events related to, for example, false negatives and false 
positives are detected by analyzing the relative relationship 
of gap fillers in the hypothesis. Thus, after determining that 

20 one of these events occurred, a cost for it can be determined 
based upon suitable statistics tables and system input 
parameters. The negative log likelihood function is desirable 
because it permits effective integration of the useful 
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components. Copies of the set of hypothetical track extensions 
which are scored are subsequently passed directly to one of the 
gate formulation and output steps 156 and 162. Note that the 
scoring process 154 also arranges the actual scores in a sparse 
5 matrix based upon observation identifiers, and passes them to 
Jc-dimensional assignment problem solving process 300. 

The assignment solving process 300 is described below. 
Its output is simply the list of assignments which constitute 
the most likely solution of the problem described by Equation 

10 [ [1]10 . [1] ]4J_. Note that both gate formulation and output 
processes 156 and 162 use (at different times) the list of 
assignments to generate the updated track history T A to 
eliminate or prune alternative previous hypotheses that are 
prohibited by the actual assignments in the list, and^ 

15 subsequently^, to output any required data. Also note that when 
one of the gate formulation and output processes 156 and 162 
accomplish these tasks, the process will subsequently generate 
and forward the new set of gates for each remaining hypothesis 
and the processes will then be prepared to receive the next set 

20 of observations. In one embodiment, the loop described here 
will generate zones for a delayed set of observations rather 
than the subsequent set. This permits processes 156 and 162 to 
operate on even observation sets while the scoring step 154 and 
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A:-dimensional solving process operate on odd sets of 
observations, or vice versa. 

The assignment solving process 300 permits the present 
invention to operate with a window size of dimension k-1 for 
5 rk>31 some k > 3 . The upper limit on k depends only upon the 
computational power of the computer 106 and the response time 
constraints of system 100. The k-1 observation sets within the 
processing window plus the prior track history result in a k- 
dimensional Assignment Problem as described by Problem 
10 Formulation [ [1]_{_Q * [1] ] 4J_. The present invention solves this 
generalized problem including the processes required to 
consider false positives and negatives, and also the processes 
required to consider sparse matrix problem formulations. 

15 I. A FIRST EMBODIMENT OF THE 

jr- DIMENSIONAL ASSIGNMENT SOLVER 300 

In describing a first embodiment of the ic-dimensional 

assignment solver 300, it is worthwhile to also discuss the 

process of Figure 4 which is used by the solver 300. Figure 4 

20 illustrates use of the Frieze and Yadagar process as shown in 

prior art for transforming a r 3 -dimensional 1 three-dimensional 

assignment problem into a r2-dimensional1 two-dimensional 

assignment problem and then use a hill climbing algorithm to 

solve the r 3-dimensionall three-dimensional assignment problem. 
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The solver 300 uses a Lagrangian Relaxation technique (well 
known in the art) to reduce the dimension of an original [k- 
dimensional 1 ic-dimensional assignment problem ( r k>31 k > 3 ) 
down to a r 3 -dimensional! three-dimensional problem and then use 
5 the process of Figure 4 to solve the r3-dimensional1 three- 
dimensional problem* Further note that the Lagrangian 
Relaxation technique is also utilized by the process of Figure 
4 and that in using this technique the requirement that each 
point is assigned to one and only one track is relaxed. 

10 Instead, an additional cost, which is equal to a respective 
Lagrangian Coefficient ["Jut"], is added to the cost or 
objective function [ f 1 (0. 11". Oil) (a) whenever a point is 
assigned to more than one track. This additional cost can be 
picked to weight the significance of each constraint violation 

15 differently, so this additional cost is represented as a vector 
of coefficients u which are correlated with respective 
observation points. Hill climbing will then develop a sequence 
of Lagrangian Coefficients sets designated (u 0 , . . . ,u 3/ 
u j+1 , . . ./Up) . That correspond to an optimum solution of the [2- 

20 dimensional 1 two-dimensional assignment problem. The 
assignments at this optimum solution are then used to "recover" 
the assignment solution of the r3-dimensional1 three-dimensional 
assignment problem. 
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In step 200 of Figure 4, initial values are selected for 
the u 0 coefficients. Because the Lagrangian Relaxation process 
is iterative, the initial values are not critical and are all 
initially selected as zero. In step 202, these additional 
5 costs are applied to the objective function f r 1 (0 . 1 r . 01 1 ) (a) . 
With the addition of the costs ["]u["] f the goal is still to 
assign the points which minimize the total cost. This 
transforms Equation r I" 1 (0 . 1 \ . 01 1 ) (a) , written for f k=3 1 k= 3 and 
altered to exclude mechanisms related to false positives and 

10 negatives, into objective function [ [2] J_l . 1 [] ]J_(a) . In the 
first iteration it is not necessary to consider the ["]u["] 
matrix because all ["]u[" ] values are set to zero. To relax 
the requirement that each point be assigned to one and only one 
track, the constraint Equation F M (0 . 1 r . 01 1 ) (d) is deleted, 

15 thereby permitting points from the last image to be assigned to 
more than one track. Note that while any axis can be chosen 
for relaxation, observation constraints are preferably relaxed. 
The effect of this relaxation is to create a new problem which 
must have the same solution in the first two axes but which can 

20 have a [differing] differ solution in the third axis. The 
result is constraints r T21 (l.lfl 1 ) (b-d) . 
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(a) [ Minimize:! Minimize Uzl z2 i i ~ u ^ih )^ z ± ± ± 

fTo l7o 123 J 3 12 



(b) Subject to: HzJ 2-J z i,i,i, =1 ' t 1 i i = If-sN^ 

(1.1) 



i 2 =i ™ 



X, z. si, [ ]i 2 = 1, . . . ,N 2J _ 

i 3 =l 1 2 3 

5 (d) 6(0,1} [ j []V z W3 ^ 

[ [2.1] ] 

Step 204 then generates from the r3-dimensional1 three- 
dimensional problem described by Problem Formulation 

10 FN ( 0 . 1 [ . 01 1 ) a new r2-dimensionall two-dimensional problem 
formulation which will have the same solution for the first two 
indices. As Problem Formulation [ [2] XL . 1 [ ] ]j_ has no 

constraints on the 3 rd axis, any value within a particular 3 rd 
axis can be used in a solution, but using anything other than 

15 the minimum value from any [3-rd] axis has the effect of 
increasing solution cost. Conceptually, the effect of step 204 
is to change the [3-dimensionall three-dimensional arrays in 
Problem Formulation [[2]XL-1[]]1 into f 2-dimensional 1 two- 
dimensional arrays as shown in Problem Formulation r [21 (1.2f11) 

20 and to generate the new r2-dimensional1 two-dimensional matrix 
2n ili2 defined as shown in Equation [ T21 (1. 3 M 1 ) . 
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(a) Minimize [: ] 



ij = 0 i 2 =0 



(c. . . -u[ ],) z 

V i, 



W, N, 



(b) Subject to: []£ £ z uu =1 ' [ ] ^ = 1 ' • • • ' N ^- 

(1.2) 



i 1= i i 3 =i 



2-1 



5 (d) I! 6 I 1 []V V 



2 



[[2.2] 



m. . = Min: arg minimize { c. . - u. I i = l, . . . , N.} 

l 1 l 2 Ijlj J L K. 



[2.3]] 



m. . = Minimum { c . . - u . I i = 1 , • . • , N. k (1.3) 

10 The cost or objective function for the reduced problem as 
defined by [ [2]J_i*2 [] ]l(a) , if evaluated at all possible values 
of u is a surface over the domain of [ U 3 ]_u. This surface is 
referred to as 0(u) and is non-smooth but provable convex 
(i.e., it has a single peak and several other critical 

15 characteristics which form terraces) . [Because of] Due to the 
convex characteristics of &(u) , the results from solving 
Problem Formulation [ [2 ] JJ, . 2 [ ] ] _L at any particular u 3 can be 
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used to generate a new set of coefficients u j+1 whose 
corresponding [Problem Formulation [2.2] problem solution is a] 
cost value is closer to the peak of <P(u) . The particular set 
of Lagrangian Coefficients that will generate the [2- 
5 dimensional! two-dimensional problem resulting in the maximum 
cost is designated u p . To recover the solution to the [3- 
dimensionall three-dimensional assignment problem requires 
solving the Equation [2JJJL.21 problem corresponding to u p . 

In step 206, the two [ ] -dimensional problem is solved 

10 directly using a technique known to those skilled in the art 
such as Reverse Auction for the corresponding cost and solution 
values. This is the evaluation of one point on the surface or 
for the first iteration <P(u 0 ) . 

Thus, after this first iteration, the points have been 

15 assigned based on all "u" values being arbitrarily set to zero. 
Because the " u" values have been arbitrarily assigned, it is 
unlikely that these assignments are correct and it is likely 
that further iterations are required to properly assign the 
points. Step 208 determines whether the points have been 

20 properly assigned after the first iteration by determining if 
for this set of assignments whether a different set of "u" 
values could result in a higher total cost. Thus, step 208 is 
implemented by determining the gradient of objective function 
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[ [2JH.2 [] ]l(a) with respect to u 3 . If the gradient is 
substantially non-zero (greater than a predetermined limit) 
then the assignments are not at or near enough to the peak of 
the [<3>(u) ] surf ace <P(u) (decision 210) , and the new set of 
5 Lagrangian Coefficients u j+1 [are]is. determined. 

Hill climbing Step 212 determines the u j+1 values based 
upon the u 3 values, the direction resulting from protecting the 
previous gradient into the U 3 domain, and a step size. The 
solution value of the r2-dimensional1 two-dimensional problem is 

10 the set of coefficients that minimize the r 2 -dimensional! two- 
dimensional problem and the actual cost at the minimum. Those 
coefficients augmented by the coefficients stored in % li2 permit 
the evaluation (but not the minimization) of the cost term in 
Problem Formulation [ [2] H . 1 [ ] ] J_. These two cost terms are 

15 lower and upper bounds on the actual minimized cost of the [3- 
dimensionall three-dimensional problem, and the difference 
between them in combination with the gradient is used to 
compute the step size. 

With this new set of ["]u["] values, steps 202-210 are 

20 repeated as a second iteration. Steps 212 and 202-210 are 
repeated until the gradient as a function of u determined in 
step 208 is less than the predetermined limit. This indicates 
that the u p values which locate the peak area of the *(u) 
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surface are determined and that the corresponding Problem 
Formulation [[2]H.2[]]1 has been solved. Step 214 will 
attempt to use the assignments that resulted from this 
particular r2-dimensionall two-dimensional assignment problem to 
5 recover the solution of the r 3 -dimensional 1 three-dimensional 
assignment problem as described below. If the limit was chosen 
properly so that the ["]u["] values are close enough to the 
peak, this recovery will yield the proper set of assignments 
that rigidly satisfies the constraint that each point be 

10 assigned to one and only one track. However, if the ["]u["] 
values are not close enough to the peak, then the limit value 
for decision 210 is reduced and the repetition of steps 212 and 
202-210 is continued. 

Step 214 recovers the f 3 -dimensional 1 three-dimensional 

15 assignment solution by using the assignment values determined 
on the last iteration through step 208. Consider the [2- 
dimensional 1 two-dimensional z assignment matrix to have 1 1 s in 
the locations specified by the list Lj={a if b^f^. If the [3- 
dimensional 1 three-dimensional z matrix is specified by placing 

20 l's at the location indicated by the list L 2 =(a i , b if m a . b .)l =1 
then the result is a solution of Problem Formulation 
[[2]J1.1[]]1. Let L 3 ={m aibi )*[ ml be the list formed by the third 
index. If each member of L 3 is unique then the L 2 solution 
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satisfies the third constraint so it is a solution to Problem 
Formulation ffl f0. ir.011). When this is not the case, recovery 
determines the minimal substitutions required within list L 3 so 
that it plus L x will be a feasible solution, i.e., a solution 
5 which satisfies the constraints of a problem formulation, but 
which may not optimize the objective function of the problem 
formulation. This stage of the recovery process is formulated 
as a r2-dimensional1 two-dimensional assignment problem [:] _as 
follows . Form a new cost matrix [c ifj ]*[, jml where c i/j =c a . bi/j for 

10 j=l. . and the N ± term is the total number of cost elements in 
the selected row of the r 3 -dimensional] three-dimensional cost 
matrix. Attempt to solve this r 2-dimensionall two-dimensional 
problem for the minimum using two constraints sets. If a 
feasible solution is found then the result will have the same 

15 form as list L x . Replace the first set of ind Texl ices by the 
indicated {a if b ± ) pairs taken from list L x and the result will 
be a feasible solution of Problem Formulation [ [1]J_Q. [1] ] 4J_. 
If no feasible solution to the new r2-dimensionall two- 
dimensional problem exists then further effort to locate the 

20 peak of <P{u) is required. 

1.1. Generalization to a Multi-Dimensional 

Assignment Solving Process 
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Let M be a fixed integer and assume that M is the 
dimension of the initial assignment problem to be solved. 
Thus, initially, the result of the scoring step 154 is a M- 
dimensional Cost Matrix which is structured as a sparse matrix 
5 (i.e., only a small percentage of the entries in the cost and 
assignment matrices are filled or non-zero) . Individual cost 
elements represent the likelihood that a track [TJTj as 
extended by the set of observations { rOvv|i=ll 0 yj - 1 i = 1 , . . .,M-1}, 
is not valid. Because the matrix is sparse the list of cost 

10 elements is stored as a packed list, and then for each 
dimension of the matrix, a vector of a variable length list of 
pointers to the cost elements is generated and stored . This 
organization means that for a particular observation [o y ]O iJ _for 
the j th list in the i th vector will be a list of pointers to all 

15 hypotheses in which [o y ] O i:/ _participates . This structure is 
further explained in the following section dealing with problem 
partitioning. 

The objective of the assignment solving process is to 
select from the set of all possible combinations of track 
20 extensions a subset that satisfies two criteria. First, each 
point in the subset of combinations should be assigned to one 
and only one track and therefore, included in one and only one 
combination of the subset, and second, the total of the scoring 
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sums for the combinations of the subset should be minimized. 
This yields the following M-dimensional equations where [k=M] k 
= M: 



*k ^ k — k 



5 (a) Minimize: v k (z k ) 



i i=0 i k =o *i-"-*k 



(b) Subject to: V, . . . X, z • • =1 i x = 1, . . .N x 



i 2 =0 i k =0 -l"*-k 



2~i - • ' 2~i 2*> • • • z i . i,. 

i 1= o ij.^o i jtl =o i k =o i k 



(c) 



for ij = 1, . . . Nj 
and j=2 , . . . k-1 



A N ^ k 



10 (d) >, ••-Z^ z i i =1 ik = 1,...N 



i^ 0 i k -i 11 



(e) i e {0,1} for all 



k 



[3.1] 

and where c k is the cost matrix [d- i]c ] which is a function of 
15 the distance between the observed point z k and the smoothed 
track determined by the filtering step, and is the cost 
function. This set of equations is similar to the set 
presented in Problem Formulation [ [1]J_0 . [1] ] 4Jl except that it 
includes the subscript and superscript k notation. Thus, in 
20 solving the M-dimensional Assignment Problem the invention 
reduces this problem to an [M-l dimensionall (Af-1) -dimensional 
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Assignment Problem and then to an [N-2 dimensional! ( A/-2 ) - 
dimensional Assignment Problem, etc. Further, the symbol 
ice{3,...,M} customizes Problem Formulation [ [3]JJL . [1] ] 4J_ to a 
particular relaxation level. That is, the notation is used to 
5 reference data from levels relatively removed as in c k+1 are the 
cost coefficients which existed prior to this level of relaxed 
coefficients c k . Note that actual observations are numbered 
from 1 to N ±J _ where N x is the number of observations in 
observation set i. Further note that the added [0]O 

10 observation in each set of observations is the unconstrained 
"gap filler." This element serves as a filler in substituting 
for missed detects, and a sequence of these elements including 
only one true observation represents the possibility that the 
observation is a false positive. Also note that by being 

15 unconstrained a gap filler may be used in as many hypotheses as 
required. 

While direct solution to [ [3]_(_!- [1] ] would give the 
precise assignment, the solution of A:-dimensional equations 
directly for large k is too complex and time consuming for 
20 practice. Thus, the present invention solves this problem 
indirectly. 

The following is a short description of many aspects of 
the present invention and includes some steps according to the 
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prior art. The first step in solving the problem indirectly is 
to reduce the complexity of the problem by the previously known 
and discussed Lagrangian Relaxation technique. According to 
the Lagrangian Relaxation technique, the absolute requirement 
5 that each point is assigned to one and only one track is 
relaxed such that for some one image, points can be assigned to 
more than one track. However, a penalty based on a respective 
Lagrangian Coefficient u k is added to the cost function when a 
point in the image is assigned to more than one track. The 

10 Lagrangian Relaxation technique reduces the complexity or 
"dimension" of the formulation of the assignment problem 
because constraints on one observation set are relaxed. Thus, 
the Lagrangian Relaxation is performed iteratively to 
repeatedly reduce the dimension until a r2-dimensionall two- 

15 dimensional penalized cost function problem results as in 
Problem Formulation [ [2 ] li. 1 [ ] ] 1. This [2-dimensional] two^ 
dimensional problem is solved then directly by a previously 
known technique such as Reverse Auction. The penalized cost 
function for the f 2-dimensionall two-dimensional problem defines 

20 a valley or convex shaped surface which is a function of 
various sets of { [u k | k=3 ] u k | k=3 . . . . , M} penalty values and one set 
of assignments for the points in two dimensions. That is, for 
each particular u 3 there is a corresponding r 2-dimensionall two- 
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dimensional penalized cost function problem and its solution. 
Note that the solution of the r 2-dimensionall two-dimensional 
penalized cost function problem identifies the set of 
assignments for the particular u 3 values that minimize the 
5 penalized cost function. However, these assignments are not 
likely to be optimum for any higher dimensional problem because 
they were based on an initial arbitrary set of u k values. 
Therefore, the next step is to determine the optimum 
assignments for the related r 3 -dimensional 1 three-dimensional 

10 penalized cost function problem. There exists a [2- 
dimensionall two-dimensional hill shaped function <P which is a 
graph of the minimums of all penalized cost functions at 
various sets of assignments in r 2 -dimension si two dimensions . 
For the f 3 -dimensional 1 three-dimensional problem, the [<i> 

15 1 function 0 can be defined based on the solution to the 
foregoing r2-dimensional1 two-dimensional penalized cost 
function. By using the current u 3 values and the { l"u k |k>31 u k | k 
> 3 1 values originally assigned [. Then], the gradient of the 
hill-shaped [<1> ] function & is determined, which [gradient] 

20 points toward the peak of the hill. By using the gradient and 
u 3 values previously selected for the one point on the hill 
(corresponding to the minimum of the penalized cost [<£> 
1 function 0 ) as a starting point, the u 3 values can be found 
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for which the corresponding problem will result in the peak of 

the ] function &. The solution of the corresponding [2- 

dimensionall two-dimensional problem is the proper values for 
f 2 1 two of the r 3 1 three sets, of indices in the [3- 
5 dimensionall three-dimensional problem. These solution indices 
can select a subsection of the cost array which maps to a [2- 
dimensionall two-dimensional array. The set of indices which 
minimize the r2-dimensional1 two-dimensional assignment problem 
based on that array corresponds to the proper assignment of 

10 points in the third dimension. The foregoing "recovery" 
process was known in the prior art, but it is modified here to 
adjust for the sparse matrix characteristic. The next task is 
to recover the solution of the proper u 4 values for the [4- 
dimensionall four-dimensional problem. The foregoing hill 

15 climbing process will not work again because the foregoing hill 
climbing process when required to locate the [4- 
dimensionan four-dimensional u 4 values for the peak of P 3 
requires the exact definition of the [$ 3 1 function (P 3 (as was 
available in the case of cP 2 ) or an always less than 

20 approximation of whereas the iteration can result in a 

greater than approximation of P 3 . According to the present 
invention, another r 3-dimensionall three-dimensional function W 
is defined which is a "less than approximation" the [3- 
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dimensional 4> 1 three-dimensional function 0 and which can be 
defined based on the solution to the [2-dimensional] two^ 
dimensional penalized cost function and the previously assigned 
and determined u k values. Next, the gradient of the [¥ 

5 ] function ¥ is found and hill climbing technique used to 

determine the u 4 values at the peak. Each selected u 4 set 
results in a new r 3 -dimensional 1 three-dimensional problem and 
requires the r 2-dimensionall two-dimensional hill climbing based 
upon new r2-dimensional1 two-dimensional problems. At the peak 

10 of the [3-dimensional W 1 three-dimensional function ¥, the 

solution values are a subset of the values required for the [4- 
dimensionall four-dimensional solution. Recovery processing 
extends the subset to a complete solution. This process is 
repeated iteratively until the u k values that result in a 

15 corresponding solution at the peak of the highest order [W and 
3> 1 functions ¥ and 0 are found. The final recovery process 
then results in the solution of ^-dimensional problem. 
Figure 5 illustrates process 300 in more detail. 

20 1.2. Problem Formulation 

In problem formulation step 310, all data structures for 
the subsequent steps are allocated and linked by pointers as 
required for execution efficiency. The incoming problem is 
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partitioned as described in the subsequent section. This 
partitioning has the effect of dividing the incoming problem 
into a set of independent problems and thus reducing the total 
workload. The partitioning process depends only on the actual 
5 cost matrix so the partitioning can and is performed for all 
levels of the relaxation process. 

1 . 2 . 1 . [ ] Relaxation and Recovery 

Step 320 begins the Lagrangian Relaxation process for 

10 reducing the M-dimensional problem by selecting all Lagrangian 
Coefficient u M penalty values initially equal to zero. The 
Lagrangian Coefficients associated with the M th constraint set 
are a[ N M +1 element! n (N^+l) -element vector. The reduction of 
this M-dimensional problem in step 322 to a [M-l 

15 dimensional! (M-l) -dimensional problem uses the two step process 
described above. First, a penalty based on the value of the 
respective u M coefficient is added to the cost function when a 
point is assigned to more than one track and then the resultant 
cost function is minimized. However, during the first 

20 iteration, the penalty is zero because all u M values are 
initially set to zero. Second, the requirement that no point 
from any image can be assigned to more than one track is 
relaxed for one of the images. In the extreme this would allow 



49 



NUMERI .01 USC 2 

a point from the relaxed image to be associate with every 
track. However, the effect of the previous penalty would 
probably mean that such an association would not minimize the 
cost. The effect of the two steps in combination is to remove 
5 a hard constraint while adding the penalty to the cost function 
so that it operates like a soft constraint. For step 322 this 
two [ Justep process results in the following penalized cost 
function problem with k=M: 



10 



(a) <P k (u k ) = Minimize: 0 k (u k ,z k ) 



i 1 = 0 



i^O 



'■k-l 



= E ...E (<...ir<)<...i. + E< 

£ . . . z * =i, i 2 = i,...n 1j _ 



(b) Subject to: 



(1.5) 



15 (c) 



(d) 



i,=0 i.^O 1^=0 i^O 



zf i 6 {0,1} 



t 



[ ]for i,- = I,. . - Nj 
[ ]and j=2, . . . k-lj_ 



[ ]_for all i z . 



k^— 



20 [ [3.2]] 
Because the constraint on single assignment of elements from 
the last image has been eliminated, a[ M-l dimensional] n (M-D- 
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dimensional problem can be developed by eliminating some of the 
possible assignments. As shown in Equations I" T 3 1 ( 1 . T31 1 6) , 
this is done by selecting the smallest cost element from each 
of the hfi h axis vectors of the cost matrix. Reduction in this 
5 manner yields a new, lower [ ]^order penalized cost function 
defined by Equations [ [3]J_1. [3] ] 6_L which has the same minimum 
cost as does the objective function defined by [ [3] XL- [2] 15) (a) 
above . 

The cost vectors are selected as follows. [ ] Define [an 



jt-i 

10 index array m-^ ^.! andja new cost array c ± mi by: 



c 



\...i k = Minimum jc*.. >ijt - u£ I ^ = 0,1, . . (1-6) 
for (i x , . . . , i k _ x ) * (0, . . . , 0) , 



k- 
C 0.. 



i u =0 I k ) 



[ [3.3]] 

The resulting [M-l dimensional] IM-l ) -dimensional problem 

is (where k=M) : 

^ k l jt l 

<P k (u k ) = Minimize: v k _ 1 (z*- 1 )= X, ... X, c i~\i k ,i k , 

15 Subject to: 2^> • • • 2^ z i i = 1' ii = 1 , ■ ■ ■ W w (1.7) 

i 2 =o i 4 =o 1 '" " 

• • • 2s 2-j • - • 2~t z i.. . .i k = ^-' 
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for ij = 
and j=2 



= 1 



f • • • 




5 




= 1,. . . 



z\ i e {0,1} 



for all i 



i • • • 



[3.4]] 



Assignment Problem [ [3]J_1. [1]]4Jl and Problem Formulation 
10 [ [3]11. [4] ] 7J_ differ' only in the dimension M vs. M-l, 
respectively. An optimal solution to [ [3]H. [4] ] 7J_ is also an 
optimal solution to equation [ [3]H. [2] ] 5J_. This relationship 
is the basis for an iterative sequence of reductions indicated 
by steps 320-332 through 330-332 and 200-204 in which the 
15 penalized cost function problem is reduced to a two [ ]- 
dimensional problem. As these formula will be used to describe 
the processing at all levels, the lowercase k is used except 
where specific reference to the top level is needed. In step 
206, the r2-dimensional1 two-dimensional penalized cost function 
20 is solved directly by the prior art Reverse Auction technique. 
Each execution of 206 produces two results, the set of z 2 values 
that minimize the problem and the cost that results v 2 when 
these z values are substituted into the objective function 



[[3]_L1. [4]]71(a) . 
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In step 208, according to the prior art, solution z 2 values 
are substituted into the r2-dimensional1 two-dimensional 
derivative of the [<I> 2 1 surface &<>. The result indicates how the 
value of u 3 should be adjusted so as to perform the hill 
5 climbing function. As was previously described the objective 
is to produce a sequence of u? values which ends when the u 3 
value is in the domain of the peak of the [0 ] surface <P. The 
section "Determining Effective Gradient" describes how new 
values are computed and how it is determined that the current 

10 Ui points to the peak of & 2 . When no further adjustment is 
required the flow moves to step 214 which will attempt to 
recover the [ 3 -dimensional 1 three-dimensional solution as 
previously described. When further adjustment is required then 
the flow progresses to step 212 and the new values of u k are 

15 computed. At the r2-dimensional1 two-dimensional level the 
method of the prior art could be used for the hill climbing 
procedure. However, it is not practical to use this prior art 
hill climbing technique to determine the updated Lagrangian 
Coefficients u k or the Max on the next (or any) higher order [<I> 

20 1 surface 0 because the next (or any) higher dimensional [<& 
1 function <P cannot be defined based on known information. 

Instead, the present invention defines a new function 
based on known information which is useful for hill climbing 
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from the third to the fourth and above dimensions, i.e., u k 
values which result in z values that are closer to the proper 
z values for the highest ic-dimension . This hill climbing 
process (which is different than that used in the prior art of 
5 Figure 4 for recovering only the three [ ] ^dimensional solution) 
is used iteratively at all lower dimensions k of the M- 
dimensional problem (including the r3-dimensional1 three- 
dimensional level where it replaces prior art) even when k is 
much larger than three. Figure 6 helps to explain this new 

10 hill climbing technique and illustrates the original k- 
dimensional cost function v k of Problem Formulation 
[ [3]J_1. [1] ] 4J_. However, the actual ^-dimensional cost surface 
v k defined by [ [3]XL- [1] ] 4J_(a) comprises scalar values at each 
point described by Ar-dimensional vectors and as such can not be 

15 drawn. Nevertheless, for illustration purposes only, Figure 6 
ignores the reality of ordering vectors and illustrates a 
concave function v^(z k ) to represent Equation [3] XI. [1]4_L- The [ 
v k ] surface v> is illustrated as being smooth to simplify the 
explanation although actually it can be imagined to be 

20 terraced. The goal of the assignment problem is to find the 
values of ~z k ; these values minimize the Jc-dimensional cost 
function v k . 
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For purposes of explanation, assume that in Figure 6, k=A 
(the procedure is used for all r k> 3 1 7c ^3 ) • This problem is 
reduced by two iterations of Lagrangian Relaxation to a [2- 

(k-l) i 

dimensional 1 two-dimensional penalized cost function 0j 
5 This cost function, and all other cost functions described 
belowj_ are also non-smooth and continuous but are illustrated 
in Figure 6 as smooth for explanation purposes. Solving the 

(k-l) 

0j 1 problem results in one set of z 2 assignments and the 
value of < ^ (k -i) i at the point u- 3 . A series of functions 

(k-l). (k-l). (k-l), (k-l), 

10 [<{)_. 1 f ... f 0^ x ] <pj , . . . 1 0i each generated from a 

different u 3 are shown. The particular series illustrates the 
process of locating the peak of < P (k -i )i ^ The r2-dimensional1 two- 
dimensional penalized cost functions 

(k_1) i (k-1)^ fMlj J k ' 1} i u i a A' 

[<J) x , ...,<|)^ ] <Pj r • • • 1 0i can be solved directly. 

15 Each such solution provides the information required to 
calculate the next u 3 value. Each iteration of the hill 
climbing improves the selection of u 3 values, i.e., yields 0 2 
problem whose solution is closer to those at the solution of 
the (p 3 problem. The result of solving 0^ is values that are 

20 on both ^ > (k . 1)1 and & k . Figure 6 illustrates the [<3> k ] surface^ 
which comprises the minimums of all k -dimensional penalized 

cost function [cj> k ] surf aces <t> k , i.e., if the Problem 

Formulations [ [3]_(JL. [2] ] 51 and [ [3]JJL. [4] ] H were solved at all 
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possible values of u k the function <P k {u k ) would result. The [G> k 
] surface & k is always less than the [v k ] surf ace_z k _ except at 
the peak as described in Equation [3]11. [5] 81 and its maximum 
occurs where the [v k 1 surf ace vy is minimum. Because the [3> k 
5 1 surface & y represents the minimum [s] of the [$ k ] sur£ace_<P k _, 
any point on the [3> k 1 surface <P y can mapped to the [v k ] surf ace 

v k _. The [ G> k ] function <P k provides a lower bound on the 

minimization problem described by [[3]H. [11311(a) . Let ~z k be 
the unknown solution to Problem Formulation [[3]JJL. [11111 and 
10 note that: 

**(u*U v k {z k )^v k ^(z k ) . (1.8) 

[[3.5]] 

Consequently, the [z k lvalues z k at the peak of <P k (i.e., the 
15 greatest lower bound on the cost of the relaxed problem) , can 
be substituted into the ^-dimensional penalized cost function 
to determine the proper assignments. Consequently, the present 
invention attempts to find the maximum on the [$ k ]surface^. 
However, it is not possible to use the prior art hill climbing 
20 to hill climb to the peak of & k because the definition of <P k 
requires exact knowledge of lower order [G> ] f unctions_^. As 
the solution of (ft is not the exact f the] solution, in that 
higher order [u lvalues of u are not yet optimized, its value 
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can be larger than the true peak of & k . As such it is not a 
lower bound on v k and it can not be used to recover the required 
solution . 

Instead, the present invention defines all auxiliary 
functions W k which [is] are based only on the solution to the 
r2-dimensional1 two-dimensional penalized cost function problem, 
lower order values z k and values [and ] u k [ values] determined 
previously by the reduction process. The [W k 1 function W k is 
a less than approximation of & kf and its gradient is used for 
hill climbing to its peak. The [z]_[ k lvalues z k at the peak of 

the [W k ] function W k are then substituted into Problem 

Formulation [ [3]JJL. [2]]5_L to determine the proper assignments. 
To define the [^ k 1 function ¥ k . the present invention explicitly 
makes the function <P k (u k ) a function of all higher order sets 
of Lagrangian Coefficients with the expanded notation: <P k (u k ; 

u k+1 , . . . , u k ) . Then, a new set of ] functions W, is defined 

recursively, using the ^'s domain [,]j_ 

¥ 3 (u 3 ) = v 2 +£ ul =<D 3 (u 3 ;u 4 , . . ., u K ) , (1.9) 

3 i 3 =o 3 

[[3.6] (a)] 
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where v 2 is the solution value for the most recent [2- 
dimensionall two-dimensional penalized cost function problem. 
For rk>31 jr >3 



,u* _1 ;u*) = 1 



,u K ), 
if known, 



f M ( U 3 ,...,u M ; U M ) + L u*, (1.10) 



otherwise. 



[ [3.6](b) 
5 ] From the definitions of <P k and v k (Problem Formulation 
t [3]H. [4] ]71 compared with Problem Formulation [ [3]H. [2] ]5JJ : 



M\ _ 



it follows that: 



W 3 (u 3 ) = v 2 +X, u^^Cu 3 ,^ 4 , 



i 3 =0 



,u M )s v 3 3 (z 3 ) 



and with that Equation [3]_Li. [5]8J_ is extended to: 



]W. (u 3 , . . . , u* -1 ; u*) <; <*V (u*;u* +1 , . . . , u M ) £ v. ( u k ) <. v. ( u *) (1. 11 



10 [ 



[3.7] ] 
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This relationship means that either $ k or W k may be used in hill 
climbing to update the Lagrangian Coefficients u. <P k is the 
preferred choice, however it is only available when the 
solution to Problem Formulation [ [3] XL. [2] ] 5J_ is a feasible 
5 solution to Problem Formulation [ [3JXL- [1] ]i_L (as in hill 
climbing from the second to third dimension which is why prior 
art worked) . For simplicity in implementation, the function W k 
is defined so that it equals <P k when either function could be 
used. It is therefore always used, even for hill climbing from 
10 the second dimension. The use of the [W 1 function ¥ which is 
based on previously assigned or determined u values and not 
higher order u values which are not yet determined, is an 
important feature of the present invention. 



15 1.2.2. Determining Effective Gradient 

After the [W k ] function iff, is defined, the next steps of 
hill climbing/peak detection are to determine the gradient of 

the [W k ] function W kf determine an increasing portion of the 

gradient and then move up the [W k ] surface W k in the direction 

20 of this increasing portion of the gradient. As shown in Figure 
6, any upward step on the [¥ k 1 surf ace W yt for example^, to the 
minimum of the [0 kj ]^ will yield a new set of values u k T 
values] (to the "left") that is closer to the ideal set of 
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values u k r values] which correspond to the minimum of the [v k 
1 function v^. While it is possible to iteratively step up this 
[¥ k 1 surface ¥ k with steps of fixed size and then determine if 
the peak has been reached, the present invention optimizes this 
5 process by determining the single step size from the starting 
point at the minimum of [^>^]^L i that will jump to the peak and 
then [calculating! calculate the values u k [ values] at the peak. 
Once the [f."u" si values u at the peak [at 1 of W k are 
determined, then the values u k [ values] can be substituted into 

10 Problem Formulation [ [3]J_1. [2] ] 5J_ to determine the proper 
assignment. (However, in a more realistic example, where k is 
much greater than three, then the values u k [ values] at the peak 
of the [W 1 function ¥ along with the lower [ 1 -order values u k [ 
values] and those assigned and yielded by the reduction steps 

15 are used to define the next higher level [¥ 1 function ¥. This 
definition of a higher order [¥ 1 function ¥ and hill climbing 
process are repeated iteratively until ¥ kr [ and] the peak of W kmLm 
and the values u k [ values] at the peak of ¥ k are identified.) 
The following is a description of how to determine the gradient 

20 of each [W ] surface ¥ and how to determine the single step 
size to jump to the peak from the starting point on each [W 
] surface ¥. 
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As noted above, each ] surf ace ¥ is non-smooth. 

Therefore, if a gradient fwerel was taken at a single point, the 
gradient may not point toward the peak. Therefore, several 
gradients (a "bundle") are determined at several points on the 

5 [ } ¥ k ] surf ace W k in the region of the starting point (i.e., 

minimum of <p k ) and then averaged. Statistically, the averaged 

result should point toward the peak of the [W k ] surf ace W k . 

Wolfe's Conjugate Subgradient Algorithm ( [Wolfe75, Wolfe79]E\. 
Wolfe. A method of conjugate subgradients for minimizing non- 

10 dif f erentiable functions. Mathematical Programming Study, 
3:145-173, 1975; P. Wolfe. Finding the nearest point in a 
polytope. Mathematical Programming Study, 11:128-149, 1976 ) for 
minimization was previously known in another environment to 
determine a gradient of a non-smooth surface using multiple 

15 subgradients and can be used with modification in the present 
invention. [Wolfe's algorithm is further described in "A 
Method of Conjugate Subgradients for Minimizing 
Nondiff erentiable Functions" page 147-173 published by 
Mathematical Programming Study 3 in 1975 and "Finding the 

20 Nearest Point in a Polytope" page 128-149 published by 
Mathematical Programming Study 11 in 1976. ] The modification 
to Wolfe's algorithm uses the information generated for 
W k (u k3 , . . . , u kk " 2 ; u kk_1 ) as the basis for calculating the 
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subgradients . 



The definition of a subgradient v of 



W k (u k3 , . . . , u kk ) is any member of the subdif f erential set defined 
as : 

5 [5W k (u) = {v6R N * +1 | (V k (u k3 , . . . , u^-^u') -W k (u k3 f . . . , u^ 1 ;^) ) >v T (u'- 
u k ) V u'ER N * +1 }] 

<?P L (u)=lv e R N fc +1 | (SUu k3 u kk ~^uM-R(u k3 u^uM ) 

>v r (u'-u k ) V u'6 R N * +1 I, 

[(]where v T is the transpose of v[)]j_ 

10 Next, a subgradient vector is determined from this 

function. If z k is the solution of Problem Formulation 
[ [3]il. [2] ] 51, then differentiating W k (u 3 , . . . , u k_1 ; u k ) with 
respect to u k k and evaluating the result with respect to the 
current selection matrix z k yields a subgradient vector: 



15 The iterative nature of the solution process at each dimension 
yields a set of such subgradients. Except for a situation 
described beloWj_ where the resultant averaged gradient does not 
point toward the peak, the most recent set of such subgradients 
are saved and used as the "bundle" for the peak finding process 
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for this dimension. For example, at the [k] level k there is 
a bundle of subgradients of the [^ k 1 surface W k near the 

minimum of the [<j) ki ] surf ace determined as a result of 

solving Problem Formulation [ [3] J_l. [1] ]4_L at all lower levels. 
5 This bundle can be averaged to approximate the gradient. 
Alternately, the previous bundle can be discarded so as to use 
the new value to initiate a new bundle. This choice provides 
a way to adjust the process to differing classes of problems, 
i.e., when data is being derived from two sensors and the 
10 relaxation proceeds from data derived from one sensor to the 
other then the prior relaxation data for the first sensor could 
be detrimental to performance on the second sensors data. 

1.2.3. [ 1 Step Size and Termination Criterion 

15 After the average gradient of the [W k ] surf ace W k is 

determined, the next step is to determine a single step that 

will jump to the peak of the [W k ] surf ace ¥ k . The basic 

strategy is [to 1 first to specify an arbitrary step size, and 
then calculate the value of W k at this step size in the 
20 direction of the gradient. If the [W k lvalue of W k is larger 
than the previous one, this probably means that the step has 
not yet reached the peak. Consequently, the step size is 
doubled and a new [W k lvalue of *F k is determined and compared 
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to the previous one. This process is repeated until the new 
[^P k lvalue of W k is less than the previous one. At that time, 
the step has gone too far and crossed over the peak. 
Consequently, the last doubling is rolled-back and that step 
5 size is used for the iteration. If the initial estimated W k is 
less than previous value then the step size is decreased by 
50%. If this still results in a smaller [W k lvalue of ¥ k , then 
the last step is rolled back and the previous step size is 
decreased by 25%. The following is a more detailed description 

10 of this process. 

With a suitable bundle of subgradients determined as just 
described [. }j_ Wolfe's algorithm can be used to determine the 
effective subgradient [(]d[)] and the [U]upgraded value Uj +1 . 
From the previous iteration, or from an initial condition, 

15 there exists a step length value [(]t[)]. The value, 

u + = Uj + td 

is calculated as an estimate of Uj +1 . To determine if the 
current step size is valid [the] we evaluate W k (u 3 , . . . , u k " 2 ; u+) . 
If the result represents an improvement then we double the step 
20 size. Otherwise we halve the step size. In either case a new 
u + is calculated. The doubling or halving continues until the 
step becomes too large to improve the result, or until it 
becomes small enough to not degrade the result. The resulting 
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suitable step size is saved with d as part of the subgradient 
bundle. The last acceptable u+ is assigned to u k +1 . 

Three distinct criteri[on]a are used to determine when u k 
is close enough to u k : 
5 1. The Wolfe's algorithm criterion of d=0 given that the test 

has been repeated with the bundle containing only the most 

recent subgradient . 

2. The difference between the lower bound <P k {u k ) and the upper 
bound v k {z k , u k ) being less than a preset relative 

10 threshold. (Six percent was found to be an effective 

threshold for radar tracking problems.) 

3. An iteration count being exceeded. 

The use of limits on the iteration are particularly 
effective for iterations at the level 3 through Jn-1[ levels, ]1 

15 as these iterations will be repeated so as to resolve higher 
order coefficient sets. With limited iterations the process is 
in general robust enough to improve the estimate of upper order 
Lagrangian Coefficients. By limiting the iteration counts then 
the total processing time for the algorithm becomes 

20 deterministic. That characteristic means the process can be 
effective for real time problems such as radar, where the 
results of the last scan of the environment must be processed 
prior to the next scan being received. 
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1.2.4. [ 1 Solution Recovery 

The following process determines if the u k values at what 

is believed to be the peak of the [W k ] function W k adequately 

approximate the u k values at the peak of the corresponding [<t> k 
5 1 function & k . This is done by determining if the corresponding 
penalized cost function is minimized. Thus, the u k values at 

the peak of the [W k ] function W k are first substituted into 

Problem Formulation [ [3]JJL . [2] ] 51 to determine a set of z 
assignments for k-1 points. During the foregoing reduction 
10 process, Problem Formulation [ [3] J_l . [4] ] 7_L yielded a tentative 
list of k selected z points that can be described by their 

indices as: [{(^i • • • ^k-l)} j-l^li** " " * where N o ^ s the 

number of cost elements selected into the solution. One 

possible solution of Problem Formulation [ [3] XI- [1] ] 4_L is the 

15 solution of Problem Formulation [ [3] XI . [2] ]5_L which is 

described as [{(i^ . . . i^m^. . . i k _ J} j° =1 ] {( *1 > • • ^-i ) } ^ 

with mf i as it was defined in Problem Formulation 
[ [3]Xi- [3] ] 6]_ . If this solution satisfies the [k-th]Jc th 
constraint set^_ then it is the optimal solution. 
20 However, if the solution for Problem Formulation 

[ [3]X1 • [2] ] 51 is not feasible (decision 355), then the 
following adjustment process determines if a solution exists 
which satisfies the [k-th]^ constraint while retaining the 
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assignments made in solving Problem Formulation r [31 (1. T41 17). 
To do this, a r2-dimensional1 two-dimensional cost matrix is 
defined based upon all observations from the [k-th]ic^ set which 
could be used to extend the relaxed solution. 

k for 1 = 0, . . . ,N k 

A J i = c i l ,... f i Jt . 1 ,i an d j = 0, . . . / Nq . (1-12) 

[[3.8]] 

If the resulting r2-dimensional1 two-dimensional assignment 



problem, 

Minimize[:] □ £ £ h ji w ji 



tt 

j=0 1=0 

t 

1=0 

\ 



Subject to: [1z2 w u = 1 j=l,...,N 0J _ 



(1.13) 



1 = 0 

[] w^e{0,l} j=0,...,N 0 , 1=0,..., N k 

[ [3.9]] 
has a feasible solution,, then the indices of that solution map 
to the solution of Problem Formulation [[3]_Q.. [1]]4J_ for the k- 
dimensional problem. The first index in each resultant 
solution entry is the pointer back to an element of the 



±1 . . . il- 1 m i ^ ^ J list. That element supplies the first k-1 
indices of the solution. The second index of the solution to 
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the recovery problem is the k th index of the solution. Together 
these indexes specify the values of z k that solve Problem 
Formulation [ [3]_[1. [1] ] 41 at the k th level. 

If Problem Formulation r [3. 91 1 (1 . 13) does not have a 
5 feasible solution^ then the value of u* which was thought to 
represent u£ is not representative of the actual peak and 
further iteration at the k th level is required. This decision 
represent the branch path from step 214 and equivalent steps. 

10 

I . 3 . Partitioning 

A partitioning process is used to divide the cost matrix 
that results from the Scoring step 154 into as many independent 
problems as possible prior to beginning the relaxation 

15 solution. The partitioning process is included with the 
Problem Formulation Step 310. The result of partitioning is a 
set of problems to be solved, i.e., there will be p 2 problems 
that consist of a single hypothesis, p 2 problems that consist 
of two hypothesis, etc. Each such problem is a group in that 

20 one or more observations or tracks are shared between members 
of the group. 

In partitioning to groups no consideration is given to the 
actual cost values. The analysis depends strictly on the basis 
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of two or more cost elements sharing the same specific axis of 
the cost matrix. In a r2-dimensional1 two-dimensional case two 
cost elements must be in the same group if they share a row or 
a column. If the two elements are in the same row, then each 
5 other cost element that is also in the same row, as well as any- 
cost elements that are in columns occupied by members of the 
row must be included in the group. The process continues 
recursively. In literature it is referred to as "Constructing 
a Spanning Forest." The £:-dimensional case is analogous to the 

10 r2-dimensionall two-dimensional case. The specific method we 
have incorporated is a depth first search, presented by Aho, 
Hopecraft, and Ullman (A.V. Aho , [in "1 J.E. Hopcroft, and J.D. 
Oilman. Design and Analysis of Computer Algorithms [, " section 
5.2, published by J^_Addison-Wes [t ] ley, MA , 1974) . 

15 The result of partitioning at level M is the set of 

problems described as ( rP ij |i=11 P i:j \ i=l , . . . , p j and j=l,...,N}, 
where N is the total number of hypothesis. The problems 
labeled { r P ±1 1 ±=1 1 P n \i=l . . . . ,p 2 } are the cases where [their] there 
is only one choice for the next observation at each scan and 

20 that observation could be used for no other track, i.e., it is 
a single isolated track. The hypothesis must be included in 
the solution set and no further processing is required. 
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As hypothesis are constructed the first element is used to 
refer to the track [id] ID. Any problem in the partitioned set 
which does not have shared costs in the first scan represent a 
case where a track could be extended in several ways but none 
5 of the ways share an observation with any other track. The 
required solution hypothesis for this case is the particular 
hypothesis with the maximum likelihood. For this case all 
likelihoods are determined as was described in Scoring and the 
maximum is selected. 

10 In addition to partitioning at the [M ] level M, 

partitioning is applied to each subsequent level M-l,...,2. 
For each problem that was not eliminated by [either by ] the 
prior selection^ partitioning is repeated, ignoring 
observations that are shared in the last set of observations.. 

15 Partitioning recognizes that relaxation will eliminate the last 
constraint set and thus partitioning is feasible for the lower 
level problems that will result from relaxation. This process 
is repeated for all levels down to k=3 . The full set of 
partitionings can be performed in the Problem Formulation Step 

20 310, prior to initiating the actual relaxation steps. The 
actual [2-dimensionall two-dimensional solver used in step 206 
includes an equivalent process so no advantage would be gained 
by partitioning at the k=2 level. 
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There are two possible solution methods for the remaining 
problems. "Branch and Bound" as was previously described, or 
the relaxation method that this invention describes. If any of 
the partitions have 5-10 possible tracks and less than 50 to 20 
5 hypotheses, then the prior art "Branch and Bound" algorithm 
generally executes faster than does the relaxation due to its 
reduced level of startup overhead. The "Branch and Bound" 
algorithm is executed against all remaining M level problems 
that satisfy the size constraint. For the remainder the 

10 Relaxation algorithm is used. The scheduling done in Problem 
Formulation allows each Problem Formulation [ [3]J_1. [2] ] 5J_ cost 
matrix resulting from the first step of relaxation to be 
partitioned. The resulting partitions can be solved by any of 
the four methods isolated track direct inclusion, isolated 

15 track tree evaluation, small group "Branch and Bound" or an 
additional stage of relaxation as has been fully described. 

The partitioning after each level of Lagrangian Relaxation 
is effective because when a problem is relaxed, the constraint 
that requires each point to be assigned to only one track is 

20 eliminated (for one image at a time) . Therefore, two tracks 
previously linked by contending for one point will be unlinked 
by the relaxation which permits the point to be assigned to 
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both tracks. The fruitfulness of partitioning increases for 
lower and lower levels of relaxation. 

The following is a more detailed description of the 
partitioning method. Its application at all but the [M ] level 
5 M depends upon the relaxation process described in this 
invention. The recursive partitioning is therefore a distinct 
part of this invention. The advantage of this method is 
greatly enhanced by the sparse cost matrix resulting from 
tracking problems. However the sparse nature of the problem 

10 requires special storage and search techniques. 

A hypothesis is the combination of the cost element c n the 
selection variable z n and all observations that made up the 
potential track extension. It can be written as, 

h n ={ c nf z n , { [ o nk =o ki | k= 1 ] Qnj^Qkj 1 k=l , . . . , [Mi] M t € { 1 , . . . , N k } } } , i.e. , 

15 cost, selection variable and observations in the hypothesis. 
rnel Here n 6 (1, . . . , AH , where N is the total number of 
hypotheses in the problem. While the cost and assignment 
matrices were previously referenced, these matrices are not 
effective storage mechanisms for tracking applications. 

20 Instead the list of all hypothesis and sets of lists for each 
dimension that reference the hypothesis set are stored. The 
hypothetical set in list form is: 
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For each dimension fk=11 k=l , . . . ,M there exists a set of lists, 
with each list element being a pointer to a particular 
5 hypothesis: 

i=N k ,k=M 

where N k is a number of hypothesis containing the [i-th]i, th 

observation from scan k. This structure is a multiply linked 
list in that any observation is associated with a set of 
pointer to all hypothesis it participates in, and any 
10 hypothesis has a set of pointers to all observations that 
formed it. (These pointers can be implemented as true pointers 
or indices depending upon the particular programming language 
utilized. ) 

Given this description of the matrix storage technique 
15 then the partitioning technique is as follows: Mark the first 
list L k . and follow out all pointers in that list to the 
indicated hypothesis h Pk . for i=l , . . . , [N ki ] . Mark all 
located hypothesis, and for each follow pointers back the 
particular [L ok ]L k for . . .^M. Those L' s if not previously 

20 marked get marked and also followed out to hypothesis elements 
and back to L's. When all such i ? s or h 1 s being located are 
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marked, then an isolated sub-problem has been identified. All 
marked elements can be removed from the lists and stored as a 
unique problem to be solved. The partitioning problem then 
continues by again starting at the first residual [L ]set_L. 
5 When none remain, the original problem is completely 
partitioned. 

Isolated problems can result from one source track having 
multiple possible extensions or from a set of source tracks 
contending for some observations. Because one of the indices 
10 of k (in our. implementation it is k=l) indicates the source 
track then it is possible to categorize the two problem types 
by observing if the isolated problem includes more than one 
L-list from the index level associated with tracks. 

15 1.4. Subsequent Track Prediction 

As noted above, after the points are assigned to the 
respective tracks, some action is taken such as controlling 
take-offs and landings to avoid collision, advising airborne 
aircraft to change course, warning airborne aircraft of an 

20 adjacent aircraft in a commercial environment, or aiming and 
firing an anti-aircraft (or anti-missile) missile, rocket or 
projectile, or taking evasive action in a military environment. 
Also, the tracking can be used to position a robot to work on 
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an object. For some or all of these applications , usually the 
tracks which have just been identified are extended or 
extrapolated to predict a subsequent position of the aircraft 
or other object. The extrapolation can be done in a variety of 
5 prior art ways; for example, straight line extensions of the 
most likely track extension hypothesis, parametric quadratic 
extension of the x versus time, y versus time, etc. functions 
that are the result of the filtering process described earlier, 
least square path fitting or Kalman filtering of just the 

10 selected hypothesis ♦ 

The process of gating, filtering, and gate generation as 
they were previously described require that a curve be fit 
through observations such that fit of observations to the 
likely path can be scored and further so that the hypothetical 

15 target future location can be calculated for the next stage of 
gating. In the implementation existing, quadratic functions 
have been fit through the measurements that occur within the 
window. A set of quadratics is used, one per sensor 
measurement. To calculate intercept locations these quadratics 

20 can be converted directly into path functions. Intercept times 
are then calculated by prior methods based upon simultaneous 
solution of path equations.' 
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The use of the fitted quadratics is not as precise as more 
conventional filters like the Kalman Filter. They are however 
much faster and sufficiently accurate for the relative scores 
required for the Assignment problem. When better location 
5 prediction is required then the assignment problem is executed 
to select the solution hypothesis and based upon the 
observations in those hypothesis the more extensive Kalman 
filter is executed. The result is tremendous computation 
savings when compared with the Kalman Filter being run on all 

10 hypothetical tracks. 

Based on the foregoing, apparatus and methods have been 
disclosed for tracking objects. However, numerous 

modifications and substitutions can be made without deviating 
from the scope of the present invention. For example, the 

15 foregoing ] functions ¥ can also be defined as recursive 

approximation problem in which several values of higher order 
u k values are used to eliminate the higher than approximation 
characteristic of the [$ k 1 function & y . The hill climbing of 
the [W ] function ¥ can be implemented by a high order hill 

20 climbing using the enhanced [$ k ] function & k . Although the 

result would not be as efficient it seems likely that the 
method would converge. Therefore, the invention has been 
disclosed by way of example and not limitation, and reference 
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should be made to the following claims to determine the scope 
of the present invention. 

II. AN ALTERNATIVE EMBODIMENT OF THE MULT I -DIMENSIONAL 

5 ASSIGNMENT SOLVING PROCESS 

The following description provides an alternative second 

embodiment of the multi-dimensional assignment solving process 

of section 1.1. However, in order to clearly describe the 

present alternative embodiment, further discussion is first 

10 presented regarding the data assignment problem of 

partitioning measurements into tracks and false alarms and, in 

particular, the representation of multi-dimensional assignment 

problems . 

15 II. 1. Formulation of the Assignment Problem 

The goal of this section is to explain the formulation of 
the data association problems, and more particularly multi- 
dimensional assignment problems, that govern large classes of 
problems in centralized or hybrid centralized-sensor level 

20 multisensor/multitarget tracking. The presentation is brief; 
technical details are presented for both track initiation and 
maintenance in JA.B. Poore[, 1 . Multi [-1 dimensional assignment 
formulation of data association problems arising from 
multitarget tracking and multisensor data fusion [, ]j_ 
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Computational Optimization and Applications , . 3[ (1994), pp. 
1 : 27-57 , 1994 ) for non-maneuvering targets and in [Ingemar J. 
Cox, Pierre Hansen, and Beta Julesz, eds., Partitioning Data 
Sets, 1 (A.B. Poore. Multidimensional assignments and 
5 multitarget tracking: Partitioning data sets. In P. Hansen, 
I.J, Cox, and B. Julesz, editors, D1MACS Series in Discrete 
Mathematics and Theoretical Computer Science , volume 19, pages 
169-198, Providence, R.I., 1995 . American Mathematical 
Society [, Providence, R.I., v. 19, Feb. 1995, pp. 169-198]! for 

10 maneuvering targets. Note that the present formulation can 
also be modified to include target features (e.g., size and 
type) -into the scoring process 154. 

The data assignment problems for multisensor and 
multitarget tracking considered here are generally posed as 

15 that of maximizing the posterior probability of the 
surveillance region (given the data) according to 



where Z M represents M data sets, y is a partition of indices of 
20 the data (and thus induces a partition of the data), T* is the 
finite collection of all such partitions, T is a discrete 




[ [5.1] ] 
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random element defined on r*, y° is a reference partition, and 
PI rr=y|z M 1 F=y\z M ) is the posterior probability of a partition y 
being true given the data Z H . The term partition is defined 
below. 

5 Consider M observation sets Z(k) (k=l, . . . , N) each of N k 

, k , N k ( k\ N k 

observations [ (z ^ / . - ] \z ± > , and let Z M denote the 

cumulative data set _defined by 

[ and Z M = {Z(l) f ...,Z(M) }, 

[5.2]] 

Z(k) =\zi\ Nk and Z m = {z(l) , . . . r Z(M)}, (2.2) 



10 respectively. In multisensor data fusion and multitarget 
tracking the data sets Z(k) may represent different classes of 
objects, and each data set can arise from different sensors. 
For track initiation the objects are measurements that must be 
partitioned into tracks and false alarms. In the formulation 

15 of track extensions a moving window over time of observations 
sets is used. The observation sets will be measurements which 
are: (a) assigned to existing tracks, (b) designated as false 
measurements, or (c) used for initiating tracks. However, note 
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that alternative data objects instead of observation sets may 
also be fused such as in sensor level tracking wherein each 
sensor forms tracks from its own measurements and then the 
tracks from the sensors are fused in a central location. Note 
5 that, as one skilled in the art will appreciate, both 
embodiments of the present invention may be used for this type 
of data fusion. 

The data assignment problem considered presently is 
represented as a case of set partitioning defined in the 
10 following way. First, for notational convenience in 

representing tracks, a zero index is added to each of the index 
sets in [ [5] ±2 . 2 [ ] ] J_, a gap filler z k 0 is added to each of the 
data sets Z(k) in [ [ 5 ] ±2 . 2 [ ] ] 1, and a "track of data" is 



15 values of 0[ and z, respectively]. A partition of the data 
will refer to a collection of tracks of data or track 
extensions wherein each observation occurs exactly once in one 
of the tracks of data and such that all data is used up. Note 
that the occurrence of the gap filler is unrestricted. The 

20 gap filler z% serves several purposes in the representation of 
missing data, false observations, initiating tracks, and 
terminating tracks. Note that a "reference partition" may be 




[and z k ] can now assume the 
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defined which is a partition in which all observations are 
declared to be false. 

Next under appropriate independence assumptions it can be 
shown that: 

5 



P(r=y°\Z M ) Y (^...i^Y 



(2.3) 



[ [5.3]] 
where L\ ± is a likelihood ratio containing probabilities for 
10 detection, maneuvers, and termination as well as probability 
density functions for measurement errors, track initiation and 
termination. Then if cj lt _ [iN= ,. ln]lN = -In (L^. . , 



-In 

15 



P(T = y\Z M ) 
p(r=y°\z M 



(i 1 ,... / i Af )6 Y N 



[[5.4] 

Using [5.3] and the zero-one variable Zi... iN = l if (iir ^ 
Y and 0] 

Using (2.3) and the zero-one variable 

20 z\ j if ( i j . , . . , ijf) 6 y, and 0, otherwise, one can then 

write the problem [ [ 5 ] J_2 . 1 [ ] ] J_ as the following M-dimensional 
assignment problem (as also presented in section Problem 
Formulation [ [3]11. [1] ]4J_ with k=M) : 
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(a) Minimize 



(b) Subject to. 



1 l = 0 i M =0 

[] S • • • £ Z i ...i =1 ' = If ■ ' ■ f^lf 



i 2 =o i M =o 



(c) [1L...E E ••■E z i i =1 - (2.5) 

for = 1 , . . . ,N k and k=2,...,M-l, 
5 (d) [] J! • • • z2 z i ...i =1 ' i M = I,- - - fN t 



(e) t]z ix...i|, 6 for ail ii/ 



AO- 



[ [5.5]] 
where c 0 ... 0 is arbitrarily defined to be zero. Here, each 

10 group of sums in the constraints f r 51 (2 . 5 f 1 1 ) (b) - 
[ [5]_L2 . 5 [] ]J_(e) represents the fact that each non-gap filler 
observation occurs exactly once in a "track of data*" One can 
modify this formulation to include multi-assignments of one, 
some, or all the actual observations. The assignment problem 

15 [ [5]J_2 • 5 [ ] ]1 is changed accordingly. For example, if z\ k is to 
be assigned no more than, exactly, or no less than n k ik times, 
then the " = 1" [in the appropriate one of the constraints 
[5.5] (b)-[5.5] (d) ]is changed to "^, =, > n k ik , " respectively. 
Expressions for the likelihood ratios L\ 1 _ Am [ such as 

20 [5.5] 1 appearing in the costs c. ^ ± =-ln IL ± ) are well [ ]- 

1 * " * M 1 " * " M 
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known. In particular, discussions of these ratios may be found 

in[: ] (_A.B. Poore[, ]_. Multi [-] dimensional assignment 

formulation of data association problems arising from 
multitarget tracking and multisensor data fusion, Computational 
5 Optimization and Applications [ , 3 (1994), pp. 27-57; and 
Ingemar J. Cox, Pierre Hansen, and Beta Julesz, eds . , 
Partitioning Data Sets, 7 . 3:27-57, 1994; A.B. Poore. 
Multidimensional assignments and multitarget tracking: 
Partitioning data sets. In P. Hansen, I.J. Cox, B. Julesz, 

10 editors, DIMACS Series in Discrete Mathematics and Theoretical 
Computer Science , volume 19, pages 169-198, Providence, R.I., 
1995 . American Mathematical Society [, Providence, R]J_. [I., v. 
19, 1995, pp. 169-198. ] Furthermore, the likelihood ratios are 
easily modified to include target features and to account for 

15 different sensor types. Also note that in track initiation, 
the M observation sets provide observations from M sensors, 
possibly all the same. Additionally note that for track 
maintenance, a sliding window of M observation sets and one 
data set containing established tracks may be used. In this 

20 latter case, the formulation is the same as above except that 
the dimension of the assignment problem is now M+l. 

II. 2. Overview of the New Lagrangian 
Relaxation Algorithms 
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Having formulated an M-dimensional assignment problem 
[ [5]J_2 . 5 [] ] we now turn to a description of the Lagrangian 
relaxation algorithm. Subsection II. 2.1 below presents many of 
the relaxation properties associated with the relaxation of an 
5 n-dimensional assignment problem to an m-dimensional one via a 
Lagrangian relaxation of n-m sets of constraints wherein [m<]m 
< n < M and preferably in the present invention embodiment 
fn-m>11 ii-2n >1 . Although any n-m constraint sets can be 
relaxed, the description here is based on relaxing the last n-m 

10 sets of constraints and keeping the first m sets. Given either 
an optimal or suboptimal solution of the relaxed problem, a 
technique for recovering a feasible solution of the 
ii-dimensional problem is presented in subsection II. 2. 2 below 
and an overview of the Lagrangian relaxation algorithm is given 

15 in subsection II. 2. 3. 

The following notation will be used throughout the 
remainder of the work. Let M be an integer such that rM>3l M > 3 
and let nE{3,...,M}. The n-dimensional assignment problem is 
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i =0 i =0 1 n 1 
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(b) Subject to[:] £ . . . ^. . .i n = 1/ ^ = 1, (2.6) 



i 2 =0 i„=0 



c /~\ 5^ ^ Z i 1 ...i n 1 

5 (c) V° wowo i n =o 1 n 

for i jc = l^ • - . t N k and k-2, . . . ,n-l, 

(d) Zi "--- i n =1 ' i n = 1 '*-*' W n' 



(e) z^.^efo,!} for all i lf ...,i fl . 

10 [6.1]] 



To ensure that a feasible solution of [ fl (2. 61". 111) always 
exists, all variables with exactly one nonzero index (i.e., 
variables of the form z n 0 . . . 0ik0 . . . 0 for i k *0) are assumed free to 
15 be assigned and the corresponding cost coefficients are well- 
defined. 



_1 1 . 2 . 1 . The Laaranqian Relaxed Assignment Problem 
The n-dimensional assignment problem [ fl (2. 6T.111) has n 
20 sets of constraints. A (N k 4-1 ) -dimensional multiplier vector 
associated with the [k-th]Jc^ constraint set will be denoted by 
u*=( u k 0f u\, . . . , u% k ) T with u k 0 =0 and k=l,...,n. The n-dimensional 
assignment problem f n (2 . 6 F . 11 1 ) is relaxed to an zn-dimensional 

85 



NUMERI .07 USC 2 



assignment problem by incorporating n-m of the n sets of 
constraints into the objective function f T 1 (2 . 6 I" . 11 1 ) (a) . 
Although any constraint set can be relaxed, the description of 
the relaxation procedure for f fl (2 . 6 I" . 11 1 ) will be based on the 
5 relaxation of the last n-m sets of constraints. The relaxed 
problem is 



[6.2]] 




u n ) = Minimize $ {z n ;u 



m+l 



(2.1) 



f • • • / 





i 1= 0 i„=0 




Subject To 





for i k = l 



AL and k = 2 



f • • • / 



^"...i/K 1 } for all i 




NUMERI .01USC2 

wherein we have added the multiplier u k 0 =0 for notational 
convenience. Thus, the multiplier u*eR Wjc+1 with u k 0 =0 is fixed. 

An optimal (or suboptimal) solution of [ [ 6] _j_2 . [2 ] ] 21 can 
be constructed from that of an m-dimensional assignment 
problem. To show this, define for each (i lf ...,i m ) an index 

(jm+l'-'-'jn) = ( Jm+lUlf • • ■ r *m) ' • • -f JnUlf • • - / ^ J ) and a neW COSt 

function d±. < by 

[ [6.3]] 

= argmin \c? ± + z2 u i I i t = 0 f . . . ,N. and k = m+l, . . . ,n\ , (2.8) 
c" ., , + £ for (i , . . .,i ) * (0, . . .0) , 



'0 . . . o ~ Z~i • • • 2-j 



Minimum 



jo,c 0 n . 



• 0i 



k=m+l 



[If ( j n ) 1 If (i ^ w . . . ,i J is not unique, choose the 

smallest such j m+lf amongst those (n-m) -tuples with the same j m+1 
choose the smallest j m +2/> e tc, so that ( j m+1 , . . . , j n ) is uniquely 
defined. ) Using the cost coefficients defined in this way, the 
following m-dimensional assignment problem is obtained: 
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gnniu** 1 , . . . ,u n ) = Minimize <p mn ( z m ;u m ^ , . . . , u n ) = vjz m ) 



= > > m m 



i, = 0 i =0 

1 m 



Subject to: - • . Si z i ...i = 1' i x = 1 , . . . ,N X , (2.9) 

A ' f\ 1 " " * ID 



i,=0 i„=0 

2 ™ 



i 1= 0 i k _ x = 0 i, +1 = 0 i m =0 

for = l,...,N k and k=2,...,m-l, 

X v • • • X j Z± , ...i — If -ijn /■ • • • 
. . - 1 in 



^...1.^(0,1} for all i lf ...,i m . 

10 [ [6.4]] 

As an aside, observe that any feasible solution z n of 

[ [112^6 [ . 1] ]1 yields a feasible solution z m of [ [6]12. [4] ] 91 via 
the construction 



m 



if Z i x ...i n Vn = 1 f ° r S ° me 

0, otherwise. 



Thus, the /n-dimensional assignment problem [[6] 12.. [4]] 91 has at 
15 least as many feasible solutions of the constraints as the 
original problem f fl f2. 6f.111) . 

The following Fact A. 2 has been shown to be true. It 
states that an optimal solution of Problem Formulation 
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[[6]J_2. [2] ]21 can be computed from that of Problem Formulation 

[ [6]_{_2 . [4] ] 91- Furthermore, the converse to this fact is 

provided in Fact A. 3. Moreover, if the solution of either of 

these two problems [ [6] 12. [2] ] H or [ [6]_L2. [4] ] 91 is e-optimal 

5 (i.e., the objective function associated with the suboptimal 

solution is within "e" of the objective function associated 

with the optimal solution), then so is the other. 

Fact A. 2. Let be a feasible solution to problem 

[ [6]12. [4] ] 91 and define w" by 

10 ^...^-w^...^ if" (i m+1 , . . .i n ) = (j m+1 , . . ■ j„) and 
(ii, ... ,ij * (0, 0) 

Wi2...i n =0 if d m+ i, ■ ■ -i n ) * (j a+ i, • • -Jn) and (2.10) 

(i 2 , . . . ,i m ) * (0, . . . , 0) 

15 

<.. 0 i ffl+I ...i n =l if cS... 0W ... in + £ u/ <0 

k=m+l 

<...oi m+1 ...i=0 if c'l.. oim ^.. in + £ u* >0^ 

k=m+l 

[[6.5]] 

20 Then w" is a feasible solution of the Lagrangian relaxed problem 
[[6]12. [2] J21 and 

0 mn {vf;u ml f . ,.,u n ) = ^d^u"* 1 , ...,u n )- £ £ u$ 

k=m+l i k =0 
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If, in addition, vf 1 is optimal for [ [ 6]_(_2. [4 ] ] 91, then vf is an 
optimal solution of [ [6]_£_2 . [2] ] 7_L and 

<p mn {u™ +i ,...,u n ) = <p; n (u- i ,...,u a )- £ i; u* 

Jt=m+1 i„=0 



With the exception of one equality being converted to an 
5 inequality, the following Fact is a converse of Fact A. 2 and 
has also been shown to be true. 

Fact A. 3. Let w" be a feasible solution to problem 
[ [6]I2 . [2] ]7_L and define w" 1 by 

iuT [] £ ■■■22 "iL.-in f ° r Ul/- • -/i-JMO, • • .,0)^. 

10 w3... 0 =0 if . . . ,i ra ) = (0, . . . ,0) and 

cS...^...^ + >0 _for all (i m+I ,. . . . (2.11) 

wS... 0 =l if (ij, . - - ,ij = (0, . . . , 0) and 

c"o...oi m+1 ...i n + £ ^0 for some ( i m+1 , . . . , i n ) . 

Jt=m+1 

[ [6.6] ] 

15 Then w* is a feasible solution of the problem [ [ 6] 12 . [ 4 ] ] 91 and 
^ n (^/u ra+1 ,...,u") > <p ma {w a tu a+1 l ... t u a )- £ u* 
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If, in addition, w" is optimal for [ [6]J_2 - [2] ] 1J_, then tv™ is an 
optimal solution of \ f 61 (2 . T4 1 1 9) , 



n+i ,...,u n )- £ i; uj 



II. 2. 2. The Recovery Procedure 
5 The next objective is to explain a recovery procedure or 

method for recovering a solution to the n-dimensional problem 
of I" M (2 . 6 1" > 11 1 ) from a relaxed problem having potentially 
substantially fewer dimensions than ["["1 (2. 61". 111). Note that 
this aspect of the alternative embodiment of the present 

10 invention is substantially different from the method disclosed 
in the first embodiment of the multi [-] dimensional assignment 
•solving process of section LI of this specification. Further, 
this alternative embodiment provides substantial benefits in 
terms of computational efficiency and accuracy over the first 

15 embodiment, as will be discussed. Thus, given a feasible 
(optimal or suboptimal) solution vf 1 of [ [6] Jj2 . [4] ] 9J_ (° r ^ if 
Problem Formulation [ [6] JJ2. [2] ]21 is constructed via Fact A. 2) , 
an explanation is provided here regarding how to generate a 
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feasible solution z n of \ \ ] (2. 61". Ill) which is close to / in a 
sense to be specified. We first assume that no variables in 
[ f 1 (2 . 6 T . 11 1 ) are preassigned to zero (this assumption will be 
removed shortly) . The difficulty with the solution is that 
5 it need not satisfy the last n-m sets of constraints in 
I'M (2. 61". Ill) . Note, however, that if W 1 is an optimal solution 
for [ [6]J_2 . [4] ] 91 and vf 1 (constructed as in Fact A. 2) satisfies 
the relaxed constraints, then is optimal for r M (2 . 6 \ . 1 1 1 ) . 
The recovery procedure described here is designed to preserve 

10 the 0-1 character of the solution vf 1 of [ [6]J_2 . [4] ] 91 as far as 
possible. That is, if fl =1 and i 2 *0 for at least one 

1=1,..., m f then the corresponding feasible solution z n of 
f fl (2 . 6 r . 11 1 ) is constructed so that x =\ for some 

( i m+1 , . . . , i n ) . Note that by this reasoning, variables of the 

15 form z n 0 _ 0im+1 _ ±n can be assigned to a value of 1 in the recovery 
problem only if However, variables Zo...oi m+1 ...i n will be 

treated differently in the recovery procedure in that they can 
be assigned 0 or 1 independent of the value This 
increases the feasible set of the recovery problem, leading to 

20 a potentially better solution. 

Let [{ (i? , . . . , i^) } N ° ] I (if, . , i^) \ N ° be an enumeration of 
J- m j = 1 l J j = i 

indices of w m (or the first m indices of w" constructed as in 
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m 

Fact A. 2) such that [w.j .j = l]w m j mj = l and (i], . . . , ijj) * 

i 1 . . .i m ^...l. 

(0,...,0). Set (ij, . . . ,i° m ) = (0, . . . , 0) for j=0 and define 



5 [ for i k =0, . . M N k ; k=m+l,...,n; j-0, . . . ,N 0 . [6.7]] 
for i k =0,...,N k ; k=m+l . . . . , n; j=0,...,N 0 . 

10 Let y denote the solution of the {n-m+1 ) -dimensional assignment 
problem: 

Minimize 2J . . . $j c£ a ** ml y i± . 

Subject to J! • • • S y-fi 1=1/ . . . ,N 0 , (2.13) 

W° ' n 



2^ ••• 2^ 2^f • • • Yji 1 ...i = i / 

7 = 0 i ,=0 i.. ,=0 i„ ,=0 i =0 111+1 n 

•* m+l Jt-1 Jt+1 n 



15 Tor i k =l,...,N k and k=m+l , . . . ,n-l , 

• • • Y-fi i 1 / — ^ f ' ' ' r N m , 

j=0 i„ +1 =0 ^=0 J -I*" » 

[ [6.8] ] 

20 The recovered feasible solution z" of f M (2 . 6 I" . 11 1 ) 
corresponding to the multiplier set {u m , . . . ,u"} is then defined 
by 
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n 



1, if . . ,,i a ) = (ii, . . for some 

j = 0, . . . ,N 'and Y L . ..i=l, (2.14) 

0, otherwise. 



[ [6.9]] 

This recovery procedure is valid as long as all cost 

coefficients c n are defined and all zero-one variables in z n are 

5 free to be assigned. Note that modifications are necessary for 

sparse problems. If the cost coefficient 

n 

[c.j . j. . }c n j j is well defined and the zero- 

one variable [z.j .j. . ] cfj tj . . is free to be 

1 1' " 'Vm+r • • n *i 

assigned to zero or one, then define 

n-m+1 n . n-m+X n 

10 2 = c.j . D • . ]c i± ± =c.j , j , . as m 

[ [6]J2. [7] ] 12J_ with z 7 n ; m+1 i being free to be assigned. 
Otherwise, z^ m+1 , is preassigned to zero. To ensure that a 
feasible solution exists, we now only need ensure that the 
variables zjo.^.o are free for j=0, 1, . . . , N 0 . (Recall that all 

15 variables of the form z n 0t _ ik _ t0 are free (to be assigned) and all 

coefficients of the form eg. . >ijt are well defined for 

k=l,...,n.) This is accomplished as follows. If the cost 
n 

coefficient [c.j.j .j n n ]c n jtj tj is well defined and 
n 1 1 1 2 * * * 1 m * * * li i 2 ■ ■ o. . .0 

[Z iW...ii 0...0 ]z ?U...iio...o is free ' then define 
12 m n 

2 « t-?0 m : 0 = c ±3i3 o...0]^"* 1 o = <V... 1 .'0...» Uith ^"'» b8in9 

1 z m 
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free. Otherwise, since all variables of the form Zo. r .± k ,, m0 are 
known feasible and have well-defined costs, put 

f n-m+1 ^ n . 1 n-m+1 ^ n 

LC jO...O " ^ c 0. . .Oijo. . .0 Jc ^—° ^ c o...oi Jk J o...o' 

5 II. 3. The Multi-Dimensional Laaranaian Relaxation 

Algorithm for the n-Dimensional Assignment 
Problem 

Starting with the M-dimensional assignment problem 
f f ] (2 . 6 r . 11 1 ) , i.e., n=M, the algorithm described below is 

10 recursive in that the M-dimensional assignment problem is 
relaxed to an m-dimensional problem by incorporating {n-m) sets 
of constrains into the objective function using Lagrangian 
relaxation of this set. This problem is maximized with respect 
to the Lagrange multipliers, and a good suboptimal solution to 

15 the original problem is recovered using an (n-m+1 ) -dimensional 
assignment problem. Each of these two (the m-dimensional and 
the {n-m+1) -dimensional assignment problems) can be solved in 
a similar manner. 

More precisely, reference is made to Figure 8 which 

20 presents a flowchart of a procedure embodying the multi- 
dimensional relaxation algorithm, referred to immediately 
above. This procedure, denoted MULTI_DIM_RELAX in Figure 8, 
has three formal parameters, n, PROB_FORMULAT ION and 
ASSIGNMT SOLUTION, which are used to transfer information 
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between recursive instantiations of the procedure. In 
particular, the parameter, n, is an input parameter supplying 
the dimension for the multi-dimensional assignment problem (as 
in [ n (2 . 6 T . 11 1 ) ) which is to be solved (i.e., to obtain an 
5 optimal or near-optimal solution) . The parameter, 

PROB_FORMULAT ION , is also an input parameter supplying the data 
structure (s) used to represent the n-dimensional assignment 
problem to be solved. Note that PROB_FORMULATION at least 
provides access to the cost matrix, [c n ] , and the observation 

10 sets whose observations are to be assigned. Additionally, the 
parameter, ASSIGNMTJSOLUTION, is used as an output parameter 
for returning a solution of a lower dimensioned assignment 
problem to an instantiation of MULTI_DIM_RELAX which is 
processing a higher dimensioned assignment problem. 

15 A description of Figure 8 follows. Assuming 

MULTI_DIM_RELAX is initially invoked with M as the value for 
the parameter n and the PROB_FORMULATION having a data 
structure (s) representing an M-dimensional assignment problem 
as in [ [5] J_2 . 5 [ ] ]]_, in step 500 an integer m, r2<m<n1 2 <m<n is 

20 chosen such that the constraint sets fori corresponding to 
dimensions jn+1, . . . , n are to be relaxed so that an m-dimensional 
problem formulation as in [ [6]_£_2 . [2] ] 7J_ results. In step 504 
an initial approximation is provided for { u m Q +1 r . . . , u n 0 ) . 
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Subsequently, in step 508 the above initial values for 
{ Uq +1 , . . . , Uq} are used in an iterative process in determining 
{uj Jm+i ,...,IP ]n ) which maximizes { & mn ( u m+i , . . . , u n ) where 
[u k EM Mk+1 ] ifeW^ll and k=m+l, . . . , n} for a feasible solution 
5 subject to the constraints of Problem Formulation 
[ [6]12. [2] ] 71. Note that by maximizing 0 mn (u m + 1 , . . . ,u n ) , some of 
the constraints that were relaxed are being forced to be 
satisfied and in so doing information is built into a solution 
of this function for solving the input assignment problem in 

10 " PROB_FORMULATION . " Further note that a non-smooth 
optimization technique is used here and that a preferred method 
of determining a maximum in step 508 is the bundle-trust method 
of Schramm and Zowe [as described in ]J_H. Schramm and J. Zowe[, 
A version of the bundle idea for minimizing a non-smooth 

15 function: Conceptual idea, convergence analysis, numerical 
results [ f 1 . SI AM Journal on Optimization, 2[ (1992)], fppl No. 
1:121-152 , February, 1992) . This method, along with various 
other methods for determining the maximum in step 508, are 
discussed below. 

20 Also note that for rm>21 m >2 , a solution to the 

optimization problem of step 508 is NP-hard and therefore 
cannot be solved optimally. That is, there is no known 
computationally tractable method for guaranteeing an optimal 
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solution. Thus, there are two possibilities: either (a) allow 
m to be greater than 2 and use auxiliary functions similar to 
those disclosed in the first embodiment of the A:-dimensional 
assignment solver 300 in section I to compute a near-optimal 
5 solution, or (b) make m=2 wherein algorithms such as the 
forward/reverse auction algorithm of D.P. Bertsekas and D.A. 
Castanon [in their paper, 1 (D.P. Bertsekas and D.A. Castanon. 
A forward/reverse auction algorithm for asymmetric assignment 
problems [, 1 . Computational Optimization and Applications , 1[ 

10 (1992), pp.]j_ 277-298, 1992) provides an optimal solution. 

If option (a) immediately above is chosen, then the 
auxiliary functions are used as approximation functions for 
obtaining at least a near-optimal solution to the optimization 
problem of step 508 . Note that the auxiliary functions used 

15 depend on the value of m. Thus, auxiliary functions used when 
m=3 will likely be different from those for jn=4 . But in any 
case, the optimization procedure is guided by using the merit 
function, <P 2n {u 3 , . . . ,u n ) , which can be computed exactly via a [2- 
dimensionall two-dimensional assignment problem for guiding the 

20 maximization process. 

Alternatively, if option (b) above is chosen, then two 
important advantages result r : 1 , namely, the optimization 
problem of step 508 can be always solved optimally and without 
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using auxiliary approximation functions. Thus, better 
solutions to the original M-dimensional assignment problem are 
likely since there is no guarantee that when the non-smooth 
optimization techniques are applied to the auxiliary functions 
5 the techniques will yield an optimal solution to step 508. 
Furthermore, it is important to note that without auxiliary 
functions, the processing in step 508 is both conceptually 
easier to follow and more efficient. 

Subsequently, in step 512 of Fig. 8, the solution 

10 (j? ]m+1 , . . . , u n ) [and w n ] is used in determining an optimal 

solution [, ] w" 1 to Problem Formulation [ [6]_L2 . [4] ] 9J_ as 

generated according to [ [6]J_2 . [3] ] 8J_ and Fact A. 3. 

In step 516, the solution vf 1 is used in defining the cost 
matrix c n " m+1 as in r f61 (2. HI 1 12) . Subsequently, if n-m+l=2,. 

15 then the assignment problem f T61 (2 . r 8 1 1 13) may be solved 
straightforwardly using known techniques such as forward/ 
reverse auction algorithms. Following this, in step 528, the 
solution to the r2-dimensional1 two-dimensional assignment 
problem is assigned to the variable ASSIGNMT_SOLUTION and in 

20 step 532 ASSIGNMT_SOLUTION is returned to a dimension three 
level recursion of MULTI_DIM_RELAX for solving a [3- 
dimensionall three-dimensional assignment problem. 
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Alternatively, if in step 520, n-m+l>2, then in step 536 
the data structure (s) representing a problem formulation as in 
r T 61 (2 . r 8 1 1 13) is generated and assigned to the parameter, 
PROB_FORMULATION . Subsequently, in step 540 a recursive copy 
5 of MULTI_DIM_RELAX is invoked to solve the lower dimensional 
assignment problem represented by PROB_FORMULATION . Upon the 
completion of step 540, the parameter, ASSIGNMT_SOLUTION, 
contains the solution to the [n-m+1 dimensional! { n-m+1) - 
dimensional assignment problem. Thus, in step 544, the [n- 

10 m+n (n-m+1) th solution is used to solve the n-dimensional 
assignment problem as discussed regarding equations 
[ [6]_L2 . [9] ] 141. Finally, in steps 548 and 552 the solution to 
the n-dimensional assignment problem is returned to the calling 
program so that, for example, it may be used in taking one or 

15 more actions such as (a) sending a warning to aircraft or sea 
facility; (b) controlling air traffic; (c) controlling anti- 
aircraft or anti-missile equipment; (d) taking evasive action; 
or (e) surveilling an object. 

20 II. 3.1. Comments on the Various Procedures 

Provided by Figure 8 

There are many procedures described by Figure 8 . One such 

procedure is the first embodiment of the multi [-] dimensional 

25 assignment solving process of section 1.1. That is, by 
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specifying m=n-l in step 500, a single set of constraints is 
relaxed in step 508, Thus, one set of constraints is. 
incorporated [ is] into the objective function via the 

Lagrangian problem formulation, resulting in an [m=n-l 

5 dimensional! (m=n-l ) -dimensional problem. The relaxed problem 
is subsequently maximized in step 512 with respect to the 
corresponding Lagrange multipliers and then a feasible solution 
is reconstructed for the n-dimensional problem using a two- 
dimensional assignment problem. The second procedure 

10 provided by Figure 8 is a novel approach which is not suggested 
by the first embodiment of section 1.1. In fact, the second 
procedure is somewhat of a mirror image of the first embodiment 
in that n-2 sets of constraints are simultaneously relaxed, 
yielding immediately an [m=2 dimensional! (m=2 ) -dimensional 

15 problem in step 512. Thus, a feasible solution to the n- 
dimensional problem is then recovered using a recursively 
obtained solution to an [n-1 dimensional! (n-1) -dimensional 
problem via step 540. In this case, the function values and 
subgradients of & 2n ( u 3 , . . . , u n ) of step 508 can be computed 

20 optimally via a two [ ] -dimensional assignment problem. The 
significant advantage here is that there is no need for the 
merit or auxiliary functions, W kf as required in the first 
embodiment of the multi [-] dimensional assignment solving 
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process of section 1.1 above and also there is no need for the 
more general merit or auxiliary functions ¥ mn , as discussed in 
subsection II. 4. 2 below. Further, note that all function 
values and subgradients used in the non[-] smooth maximization 
5 process are computed exactly (i.e., optimally) in this second 
procedure. Moreover, problem decomposition is now carried out 
for the n-dimensional problem; however, decomposition of the 
[n-1 dimensional! (n-1) -dimensional recovery problem (and all 
lower recovery problems) is performed only after the problem is 

10 formulated. Between these two procedures are a host of 
different relaxation schemes based on relaxing n - m sets of 
constraints to an m-dimensional problem ( T2<m<n1 2<m<n ) , but 
these all have the same difficulties as the procedure for the 
first embodiment of section 1.1 in that the relaxed problem is 

15 an NP-hard problem. To resolve this difficulty, we use an 
auxiliary or merit function W mn as described in subsection 
II. 4.2 below. (The notation W k was used for ¥ n _ lfli in section 
1.1.) For the case jn<n-l, the recovery procedure is still 
based on an NP-hard ( rn-m+l-dimensionall n-m+1) -dimensional 

20 assignment problem. Note that the partitioning techniques 
similar to those discussed in Section 1.3 may be used to 
identify the assignment problem with a layered graph and then 
to find the disjoint components of this graph. In general, all 
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relaxed problems can be decomposed prior to any non[-] smooth 
computations because their structure stays fixed throughout the 
algorithm of Fig. 8. However, all recovery problems cannot be 
decomposed until they are formulated, as their structure 
5 changes as the solutions to the relaxed problems change. 

II. 4. Details and Refinements Relating to the 

Flowchart of Figure 8 

Further explanation is provided here on how various steps 

10 of Figure 8 are solved. Note that the refinements presented 

here can significantly increase the speed of the relaxation 

procedure of step 508. 

II. 4.1 Maximization of the Non-smooth Function 

One of the key steps in the Lagrangian relaxation 
algorithm in section II. 3 is the solution of the problem 

Maximize {4> mn (u m+1 , . . . ,u n ) : fu k eR Mk+1 1 u k eR"*+ j ; [k=m+l] jc= 
20 m+1 . . . . , n}j_ 

[ f7.n (2.15) 

where u k 0 =0 for all k=m+l,...,n. The evaluation of 

0 mn ( u m+1 , . . . , u n ) requires the optimal solution of the 
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corresponding minimization problem [[6]J_2. [2] Z n varies for each 
instance of (u m+1 , . . . ,u n ]l) . [ ] The following discussion provides 
some properties of these functions. 

Fact A. 4 . Let u m+1 , . . . , u n be multiplier vectors associated 
5 with the (m+l) st through the n th set of constraints 
[ []12^6[.l] ]1, let 0 mn be as defined in [ [ 6] J_2 . [2 ] ] 71, let 
Vnti 7 3n ) be the objective function value of the n-dimensional 
assignment problem in equation I" M (2 . 6 T . 1 1 1 ) , let z n be any 
feasible solution of [fl (2. 6 r . 11 1 ) , and let J 7 ]n be an optimal 
10 solution of r M (2. 6 T . Ill ) . Then, <P mn ( u m+1 , . . . , u n ) is piecewise 
affine, concave and continuous in { u m+1 , . . . , u n } and 

& ^(u m+1 u n ) < V n (I n ) * V n (z"± (2.16) 

[[7.2a]] 

15 Furthermore, 

<P m . ln {u\u m+1 ,...,u n ) < & mn {u m+1 , . . . ,u nl) 

[7.2a] 

for all uk£MMk+l with v=0 and k=m, . . . ,n. ] j (2 17) 

for all u k €R Mk+1 with u ft =0 and k=m n. 

20 

Most of the procedures for non[-] smooth optimization are 
based on generalized gradients called subgradients , given by 
the following definition. 
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Definition . At u=( u m+1 , . . . , u n ) the set [ 5<P mn l d@ mn (u) is called 

a subdifferential of & mn and[ is] defined by 

[od> mn ] 

d0 mn (u)={ rqgR^^xl q g R«*n +1 x . . . r xM Mn+1 1 3^ 1 xR^| <p_ , ( w) [-<I> ran ]- 
5 0 ran (u) [<.q T ]ZgL ( [w-u)Vw6l Mm+1+1 x. . . xM Mn+1 } 

[7.3] 

A vector qg5$ mn (u) 1 w-u) (2.18) 

W € R^Wx. . .xR"n+*). 

10 A vector q € d<t> m n (u) is called a subgradient . 

If z" is an optimal solution of [ [6] 12. . [2] ] 7_L computed 
during evaluation of <P mn (u) , differentiating <P mn with respect to 
u" n yields the following [i n -th]X," component of a subgradient 
g of <P mn (u) 

15 g 0 *=0, 

<=E---E E --E ^...i- 1 (2.19) 

i 1 = 0 ^=0 i ktl =0 i n =0 

for i k =l,...,N k and k=m+l,...,n. 

[[7.4]] 

If z n is the unique optimal solution of r f 61 (2 . T21 1 7 ) . 
20 [6<l> mn (u) ={g} , and <I> mn l d0 m n I u) = f a) and is dif f erentiable at u. 
If the optimal solution is not unique, then there are finitely 
many such solutions, say z n ( 1 ),..., z" {K) . Given the 
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corresponding subgradients, q l ,...,<f t the subdif f erential 
[53>]^£(u) is the convex hull of { g 1 , . . . , g*} . 



II. 4.2 Mathematical Description of 
5 a Merit or Auxiliary Function 

For real-time needs, one must address the fact that the 

non-smooth optimization problem of step 508 (Fig. 8) requires 

the solution of an NP-hard problem for m>2 . One approach to 

this problem is to use the following merit or auxiliary 

10 function to decide whether a function value has increased or 

decreased sufficiently in the line search or trust region 

methods : 

W m Ju 3 , . . . ,u m ;u m+1 , . . . ,u n ) = (2.20) 

®„(u" +1 ,...,0, ifm = 2 

or (2.7) is solved optimally, 
® 2n ( u 3 , . . . , u m ; u m+1 , . . . , u n ) , otherwise. 

15 

[[7.5]] 

where the multipliers jF, . . . , if that appear in lower order 
relaxations used to construct (suboptimal) solutions of the 
/n-dimensional relaxed problem ([6] 2. [2] 7.) have been explicitly 
20 included. Note that W mn is well-defined since ([6]2. [4] 9) can 
always be solved optimally if m=2 . For sufficiently small 
problems ([6]2. [2]7J or ( [6] 2. [4].9) , one can more efficiently 
solve the NP-hard problem by branch and bound. This is the 
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reason for the inclusion of the first case; otherwise, the 
relaxed function <P 2n is used to guide the non[-] smooth 
optimization phase. That the merit function provides a lower 
bound for the optimal solution follows directly from Fact A. 4 
5 and the following fact. 

Fact A. 5. Given the definition of W mn in ([7] 2. [5] 20) 
above_^ 

W mn {}P,.. . f ^;\^ +1 f ... f u n ) [<® mn ]£® mn (u m+1 ,...,u n ) for all multipliers 

10 j?, . . . ^u 5 ,^, . . . ,u n . (2.21) 

[[7.6]] 

The actual function value used in the optimization phase is W mn ; 
however, the subgradients are computed as in ( [7 ] 2. [ 4 ] 19.) , but 
with the solution z\ 1 ^_ in being a suboptimal solution constructed 

15 from a relaxation procedure applied to the jn-dimensional 
problem. Again, the use of these lower order relaxed problems 
is the reason for the inclusion of the multipliers j?,...,^. 

To explain how the merit function is used, suppose there 
is a current multiplier set ( u£J^ . . . , u n old ) and it is desirable 

20 to update to a new multiplier set ( u%l, - . . , u n new ) via 
(uKi,.../"ZeJ = (uS!^...,uSid) + (Axf* 1 ,...,*}?) . Then 
^mn (^id/ - • • r "Sid/ u Tid f - • • / u n old ) is computed where ( ]? old , . . . , u m old ) is 
obtained during the relaxation process used to compute a high 

107 



NUMERI .01 USC 2 

quality solution to the relaxed jn-dimensional assignment 
problem ([6] 2. [2] 7) at (u m 'V..,u n ) - ( uf 16 , . . . , u n old ) . The 
decision to accept ( u%l, . . . , u" ew ) is then based on 

^mn (-Moid' • • * r U old' ^oid' • • • r U old) ^ ^mn (MneW * • * r iw' ^new f * * * f U jiew) Or 

5 some other stopping criteria commonly used in line searches. 
Again, {u^ ew , . - . , uS ew ) represents the multiplier set used in the 
lower level relaxation procedure to construct a high quality 
feasible solution to the ^-dimensional relaxed problem 
([6]2.[2]7) at (u m+1 , . . . ,u n ) = (u m ne l, . . - ,u n new ) . The point is that 

10 each time one changes ( u ml , . . . , u n ) and uses the merit function 
Van (IP/ • • • / - • • / u n ) for comparison purposes, one must 

generally change the lower level multipliers {U 3 , . . . , if 1 ) . 

An illustration of this merit function for [m=n-l is given 
in "PARTITIONING MULTIPLE DATA SETS, MULTIDIMENSIONAL 

15 ASSIGNMENTS AND LAG RANG I AN RELAXATION", by Aubrey B. Poore and 
Nenad Rijavec, and in "QUADRATIC ASSIGNMENT AND RELATED 
PROBLEMS, DIMACS SERIES IN DISCRETE MATHEMATICS AND THEORETICAL 
COMPUTER SCIENCE", in American Mathematical Society, 
Providence, R.I., Vol. 16, 1994, pp. 25-37 edited by Panos 

20 1 m=n-l is given in (A.B. Poore and N. Riiavec. Partitioning 
multiple data sets: multidimensional assignments and Lagrangian 
relaxation. In P. M. Pardalos and [Henry Wolkowicz. 
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l H. Wolkowicz, editors, Oqudratic assignment and related 
problems : DIMACS Series in Discrete Mathematics and Theoretical 
Computer Science, volume 16, pages 25-37, 1994). 

5 II. 4.3 Non-smooth Optimization Methods 

By Fact A. 4 the function & mn (u) is a continuous, piecewise 
affine, and concave, so that the negative of @ mn (u) is convex. 
Thus the problem of maximizing <P mn (u) is one of non[-] smooth 
optimization. Since there is a large literature on such 

10 problems, only a brief discussion of the primary classes of 
methods for solving such problems is provided. A first class 
of methods, known as subgradient methods, are reviewed and 
analyzed in [the book by ]J_N.Z. Shor[,]_^ Minimization Methods 
for Non~Differentiable Functions[ , ]^ Springer-Verlag, New York, 

15 19851. Despite their relative simplicity, subgradient methods 
have some drawbacks that make them inappropriate for the 
tracking problem. They do not guarantee a descent direction at 
each iteration, they lack a clear stopping criterion, and 
exhibit poor convergence (less than linear) . 

20 A more recent and sophisticated class of methods are the 

bundle methods; excellent developments are presented in [the 
books of Hiriart-Urruty and Lemarechal: ]JJ.-B. Hiriart-Urruty 
and C. Lemarechal [, Convex Analysis and Minimization 
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Algorithms I and Springer-Verlag, Berlin, 1993; [and the 
book of ]K.C. Kiwiel [ : ]_. Methods of Descent for Non- 
Different iable Optimization . In A. Dold and B. Eckmann , 
editors. Lecture Notes in Mathematics, l"Vol . 1 volume 1133, 
5 Berlin, 1985. Springer-Verlag \ , Berlin, 1985]!. Bundle methods 
retain a set (or bundle) of previously computed subgradients to 
determine the best possible descent direction at the current 
iteration. Because of the non [-] smoothness of <P mnr the 
resulting direction may not provide a "sufficient" decrease in 

10 0 mn . In this case, bundle algorithms take a "null" step, 
wherein the bundle is enriched by a subgradient close to u k . 
As a result, bundle methods are non-ascent methods which 
satisfy the relaxed descent condition @ mn (u k+1 ) [<® mn ] <® m n ( u k) if 
[u k+1 *u k ] u k+1 *u k . These methods have been shown to have good 

15 convergence properties. In particular, bundle method variants 
have been proven to converge in a finite number of steps for 

piecewise affine convex functionals in _[K.C. Kiwiel [, ]_. An 

[AJaggregate [S] subgradient [M]method for [N] non-smooth 
[CJconvex [M] minimization [ , 1 . Mathematical Programming^ 27 [, 

20 (1983) pp. ]i.320-341, Fandl 1983; H. Schramm and J. Zowe[, ]^_A 
version of the bundle idea for minimizing a non-smooth 
function: Conceptual idea, convergence analysis, numerical 
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results [, 1 . SI AM Journal on Optimization, 2[ (1992)], fppT No . 
1:121-152 , February, 1992) . 

All of the above non-smooth optimization methods require 
the computation of_& mn {u) and a rqe5<t> mn l a € d<P m n (u) . These in 
5 turn require the computation of the relaxed cost coefficients 
[ [6]J_2 . [3] ] 8J_. In both the first and second procedures 
discussed in section II. 3.1, maximization of @2n( u ) must be 
repeatedly evaluated. In the most efficient implementations 
presently known of these two procedures, it was found that at 

10 least 95% of the computational effort of the entire procedure 
is spent in the evaluation of the relaxed cost coefficients 
[ [6]_£_2 . [3] ] 8J_ as part of computing ^( u )* Thus, generally a 
method should be chosen that makes as efficient use of the 
subgradients and function values as [is ] possible, even at the 

15 cost sophisticated manipulation of the subgradients. In 
evaluating three different bundle procedures: (a) the conjugate 
subgradient method of Wolfe used in section 1.1 of the first 
embodiment of the present invention; (b) the aggregate 
subgradient method of Kiwiel ([described in ]K.C. Kiwiel[, ]_,_ 

20 An [A] aggregate [S] subgradient [M] method for [N] non-smooth 
[CJconvex [M] minimization [ , 1 . Mathematical Programming^ 27 [, 
(1983) pp. 1 : 320-341 , 1983 ) ; and (c) the bundle trust method of 
Schramm and Zowe ([described in ]H. Schramm and J. Zowe[, 1 . A 
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version of the bundle idea for minimizing a non-smooth 
function: Conceptual idea, convergence analysis, numerical 
results [, 1 . SI AM Journal on Optimization, 2[ (1992)], [pp]No. 
1:121-152 , February, 1992 ) , it was determined that for a fixed 
5 number of non-smooth iterations, say, ten, the bundle-trust 
method provides good quality solutions with the fewest number 
of function and subgradient evaluations of all the methods, and 
is therefore the preferred method. 

10 II. 4. 4 The Two F 1 -Dimensional Assignment Problem 

The forward/reverse auction algorithm of Bertsekas and 
Castanon ([as described in ]D.P. Bertsekas and D.A. Castanon[, 

]_. A forward/reverse auction algorithm for asymmetric 

assignment problems [, ]_. Computational Optimization and 

15 Applications, 1[ (1992), pp. 1 :277-298 , 1992 ) is used to solve 
the many relaxed two [ ] -dimensional problems that occur in the 
course of execution. 

II. 4. 5. Initial Multipliers and Hot Starts 
20 The effective use of "hot starts" is fundamental for real- 

time applications. A good initial set of multipliers can 
significantly reduce the number of non-smooth iterations (and 
hence the number of & mn evaluations) required for a high quality 
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recovered solution. Further, there are several ways that 
multipliers can be reused. First, if an n-dimensional problem 
is relaxed to an m-dimensional problem, relaxation provides the 
multiplier set {if 1 * 1 , u m+2 , . . . , u n } . These can be used as the 
5 initial multipliers for the [n-m+1 dimensionall (n-m+D- 
dimensional recovery problem for n-jn+l>2. This approach has 
also worked well to reduce the number of non-smooth iterations 
during recovery. 

Further, for track maintenance, initial feasible solutions 

10 are generated as follows. When a new scan of information (a 
new observation set) arrives from a sensor, one can obtain an 
initial primal feasible solution by matching new reports to 
existing tracks via a two [ ] -dimensional assignment problem. 
This is known as the track-while-scan (TWS) approach. Thus, an 

15 initial primal solution exists and then we use the above 
relaxation procedure [is ]to make improvements to this TWS 
solution. Also for track maintenance, multipliers are 
available from a previously solved and closely related multi[- 
] dimensional assignment problems for all but the new 

20 observation set. 
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II. 4. 6. Local Search Methods 
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Given a feasible solution of the multi [-] dimensional 
assignment problem, one can consider local search procedures to 
improve this result, as described in J.C.H. Papadimitriou and K. 
Steiglitz [ , Combinatorial Optimization: Algorithms and 
5 Complexity[ , ] ^ Prentice-Hall , Inc. , Englewood Cliffs, NJ, 
1982[, and]j__ J. Pearl f,]^ Heuristics : [ i] Intelligent [sJSearch 
[s] Strategies for [cJComputer [p] Problem [ s] Solving[ , ] ^ 
Addison-Wesley, Reading, MA, 1984) . One method is the idea of 
[ k-interchanaesl k interchanges . Since for sparse problems only 

10 those active arcs that participate in the solution are stored, 
it is difficult to efficiently identify eligible variable 
exchange sets with some data structures for solving the 
assignment problem. However, a local search that may be more 
promising is to use the two [ ] ^dimensional assignment solver in 

15 the following way. Given a feasible solution to the multi [- 
] dimensional assignment problem, the indices that correspond to 
active arcs in the solution are enumerated. Subseguently, one 
coordinate is freed to formulate a two- dimensional assignment 
problem with one index corresponding to the enumeration and the 

20 other to the freed coordinate, and then solve a two [ ]- 
dimensional assignment problem to improve the freed index 
position . 
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II. 4. 7. Problem Decomposition 
The procedures described thus far are all based on 
relaxation. Due to the sparsity of assignment problems, 
however, frequently a decomposition of the problem into a 
5 collection of disjoint components can be done wherein each of 
the components can be solved independently. Due to the setup 
costs of Lagrangian relaxation, a branch and bound procedure is 
generally more efficient for small components, say ten to 
twenty feasible 0-1 variables (i.e., z± i) . Otherwise, the 

10 relaxation procedures described above are used. Perhaps the 
easiest way to view the decomposition method is to view the 
reports or measurements as a layered graph. A vertex is 
associated with each observation point, and an edge is allowed 
to connect two vertices only if the two observations belong to 

15 at least one feasible track of observations. Given this graph, 
the decomposition problem can then be posed as that of 
identifying the connected subcomponents of a graph which can be 
accomplished by constructing a spanning forest via a depth 
first search algorithm, as described in J.A.V. Aho, J.E. 

20 Hopcroft, and J.D. Ullman[,]^_ The Design and Analysis of 
Computer Algorithms [,] Addison-Wesley [ Publishing Company, 
Reading], MA, 19741. Additionally, decomposition of a 

different type might be based on the identification of bi- and 
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tri-connected components of a graph and enumerating on the 
connections. Here is a technical explanation. Let 

2 = {Zi 1 i 2 ...[in/ziii2]dLZi 1 i 2 ...i n is n °t preassigned to zero} 
denote the set of assignable variables. Define an undirected 
5 graph G{N,A) where the set of nodes is 

AT={ rz n |n=11 z n ± \n=l , . . . ,n; i n =l, . . . , NJ 

and arcs, 

A={ Iz^^z 1 ^) r | n^l 1 \n*l , j n *0, j x *0 and there exists z ±li2 _^e% 

such that j n =i n and ji=ij}. 

10 v Note that the nodes corresponding to zero index have not 

been included in the above defined graph, since two variables 
that have only the zero index in common can be assigned 
independently. Connected components of the graph are then 
easily found by constructing a spanning forest via a depth 

15 first search. (A detailed algorithm can be found in the book 
by Aho, Hopcroft and Ullman cited above) . Furthermore, this 
procedure is used at each level in the relaxation, i.e., is 
applied to each assignment problem [ [3] J_l. [1] ] 4J_ for n=3, . . . ,M. 
The original relaxation problem is decomposed first. All 

20 relaxed assignment problems can be decomposed a priori and all 
recovery problems can be decomposed only after they are 
formulated. Hence, in the n-to-(n-l) case, we have n-2 relaxed 
problems that can all be decomposed initially, and the recovery 
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problems that are not decomposed (since they are all two [ ]- 
dimensional) . In the n-to-2 case, we have only one relaxed 
problem that can be decomposed initially. This case yields n-3 
recovery problems , which can be decomposed only after they are 
5 formulated. 

III. NEW APPROACHES TO TRACK INITIATION AND MAINTENANCE 
USING MULTIDIMENSIONAL ASSIGNMENT PROBLEMS 

10 III . 1 . Introduction 

The ever-increasing demand for sensor surveillance systems 
is to accurate [ly] track and identify objects in real-time, 
even for dense target scenarios and in regions of high track 
contention. The use of multiple sensors, through more varied 

15 information, has the potential to greatly improve target state 
estimation and identification. However, to take full advantage 
of the available data, a multiple frame data association 
approach is needed. The most popular such approach is an 
enumerative technique called multiple hypothesis tracking 

20 (MHT) . As an enumerative technique, MHT can be too 
computationally intensive for real-time needs because this 
problem is NP-hard. A promising approach is to utilize the 
recently developed Lagrangian relaxation algorithms 
(K. P. Pattipati, S. Deb, and Y . Bar-Shalom. A multisensor- 

25 multitaraet data association algorithm for heterogeneous 
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sensors. IEEE Transactions on Aerospace and Electronic Systems, 
29, No. 2:560-568, April 1993; A.B. Poore and N.Riiavec. 
Partitioning multiple data sets: multidimensional assignments 
and lagrangian relaxation, in quadratic assignment and related 
5 problems, P.M. Pardalos and H.Wolkowicz, editors, DIMACS series 
in Discrete Mathematics and Theoretical Computer Science, 
16:25-37, 1994; A.J. Robertson III. A class of lagrangian 
relaxation algorithms for the multidimensional problem. Ph.D. 
Thesis, Colorado State University* Ft. Collins, CO, 1995) for 

10 the multidimensional assignment problem; however, there are 
many other potentially good approaches to these assignment 
problems such as LP relaxation combined with an interior point 
method, GRASP, and parallelization. 

These data association problems are fundamentally 

15 important in tracking. Thus, our first objective is to present 
an overview of the tracking problem and the context within 
which these data association problems fit. Following a brief 
review of the probabilistic framework for data association, the 
second objective is to formulate two models for track 

20 initiation and maintenance as multidimensional assignment 
problems. This formulation is based on a window moving over 
the frames of data. The first and simpler model uses the same 
window length for track initiation and maintenance, while the 
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second model uses different lengths. As the window moves over 
the frames of data, one creates a sequence of these 
multidimensional assignment problems and there is an overlap 
region of frames common to adjacent windows. Thus, given the 
5 solution of one problem, one can develop a warm start for the 
solution to the next problem in the sequence, as shown 
hereinafter. Such information is critically important to the 
design of real-time algorithms. 

[ ] 

10 III. 2. Overview of the Tracking Problem 

Tracking and data fusion are complex subjects of 
specialization that can only be briefly summarized as they are 
related to the subject of this paper. The processing of track 
multiple targets using data from one or more sensors is 

15 typically partitioned into two stages or major functional 
blocks, namely, signal processing and data r process 1 processing 
( Y. Bar-Shalom. Multitaraet-Multisensor Tracking: Advanced 
Applications. Artech House, Dedham, MA, 1990; Y. Bar-Shalom and 
T.E. Fortmann. Tracking and Data Association. Academic Press, 

20 Boston, MA, 1988; S.S. Blackman. Multiple Target Tracking with 
Radar Applications . Artech House, Dedham, MA, 1986) . The first 
stage is the signal processing that takes raw sensor data and 
outputs reports. Reports are sometimes called observations, 
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threshold exceedances, plots, hits, or returns, depending on 
the type of sensor. The true source of each report is usually 
unknown and can be due to a target of interest, a false signal, 
or persistent background objects that can be moving in the 
5 field of view of the sensor. 

Two principal functions of data processing are tracking 
and discrimination or target identification. The 
discrimination or target identification function distinguishes 
between tracks that are for targets of interest and other 

10 tracks such as those due to background. Also, if there is 
enough information, each track is classified as to the type of 
target it appears to represent. Target identification and 
discrimination will not be discussed further in this paper 
except to comment here on the use of attributes (including 

15 identification information) in the data association. 

There are two types of attributes or features, namely, 
discrete valued variables and continuous valued variables. 
Available attributes of either or both types should be used in 
computing the log likelihoods or scores for data association. 

20 In discussing tracking in this paper, the attributes will not 
be mentioned explicitly but it should be understood that some 
reports and some tracks may include attributes information that 
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should be and can be used in the track processing (both data 
association and filtering) if it is useful. 

[ ] 

III. 2.1. Tracking Functions 

5 The tracking function accepts reports for each frame of 

data and constructs or maintains a target state estimate, the 

variance-covariance matrix of the state estimation error, and 

refined estimates (or probabilities) of the target attributes. 

The state estimate typically includes estimates of target 
10 position and velocity in three dimensions and possibly also 

acceleration and other variables of interest. 

The tracking function is typically viewed in two stages, 

namely, data association and filtering. Also, the process of 

constructing new tracks, called track initiation, is different 
15 from the process of updating existing tracks, called track 

maintenance . 

In track maintenance, the data association function 
decides how to relate the reports from the current frame of 
data to the previously computed tracks. In one approach, at 
20 most one report is assigned to each track, and in other 
approaches, weights are assigned to the pairings of reports to 
a track. After the data association, the filter function 
updates each target state estimate using the one or more (with 
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weights) reports that were determined by the data association 
function. A filter commonly used for multiple target tracking 
and data fusion is the well known Kalman filter (or the 
extended Kalman filter in the case of nonlinear problems) or a 
5 simplified version of it (Y. Bar-Shalom. Multitarget- 
Multisensor Tracking: Advanced Applications . Artech House, MA, 
1990; Y. Bar-Shalom and T.E. Fortmann. Tracking and Data 
Association. Academic Press, Boston, MA, 1988; S.S. Blackman. 
Multiple Target Tracking with Radar Applications . Artech House, 

10 Dedham, MA, 1986) . 

In track initiation, typically a sequence of reports is 
selected (one from each of a few frames of data) to construct 
a new track. In track initiation, the filtering function 
constructs a target state estimate and related information 

15 based on the selected sequence of reports. The new track is 
later updated by the track maintenance processing of the 
subsequent frames of reports. 

In some trackers, there is also a tentative tracking 
function for processing recently established tracks until there 

20 is enough confidence to include them in track maintenance. 
While for simplicity of presentation, tentative tracking will 
not be included in the processing discussed in this paper, the 
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techniques discussed could readily include one or more 
tentative tracking functions. 

There have been numerous approaches developed to perform 
the data association function. Since optimal data association 
5 is far too complex to implement, good but practical sub-optimal 
approaches are pursued. Data association approaches can be 
classified in a number of ways. One way to classify data 
association approaches is based on the number of data frames 
used in the association process \ . 1 (O.E. Drummond. Multiple 

10 sensor tracking with multiple frame, probabilistic data 
association. In Signal and Data Processing of Small Targets , 

SPIE Proceedings , volume 2561, pages 322-336, 1995) . In 

single frame data association for track maintenance, "hard 
decisions" are made [each frame as to how the"! on the assignment 

15 of reports [ are to be related] to the tracks. Some single 
frame approaches include: individual nearest neighbor, PDA, 
JPDA, and global nearest neighbor (sequential most probable 
hypothesis tracking), which uses a two [ ] -dimensional 
assignment algorithm (Y. Bar-Shalom. Multitarget-Multisensor 

20 Tracking: Advanced Applications . Artech House, MA, 1990; Y. 
Bar-Shalom and T.E. Fortmann. Tracking and Data Association. 
Academic Press, Boston, MA, 1988; S.S. Blackman. Multiple 
Target Tracking with Radar Applications . Artech House, Dedham, 
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MA, 1986) . In most single frame data association approaches, 
only one track per object is carried forward to be used in 
processing the next frame of reports. 

Multiple frame data association is more complex and 
5 frequently involves purposely fcarryl carrying forward more than 
one track per target to be used in processing the next frame of 
reports. By retaining more than one track per target, the 
tracking performance is improved at the cost of increased 
processing load. The best known multiple frame data 

10 association approach is an enumerative technique called 
multiple hypothesis tracking (MHT) which enumerates all the 
global hypothesis with various pruning rules. 

As shown hereinafter, the log likelihood of the 
probability of a global hypothesis can be decomposed into the 

15 sum of scores of each track that is contained in the 
hypothesis. As a consequence, the most probable hypothesis can 
be identified with the aid of a non-enumerative algorithm for 
the assignment so formulated. 

[ ] III. 2. 2. Interface Between Track Initiation 

20 And Track Maintenance 

There are a number of ways that track initiation and track 

maintenance functions can interact (O.E . Drummond. Multiple 

target tracking with multiple frame, probabilistic data 

association. In Signal and Data Proceedings , volume 1954, pages 
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394-408, 1993) . Two ways that are especially pertinent to the 
methods of this paper will be discussed in this section. 

One approach is to treat these two functions sequentially, 
by first assigning reports to tracks and then using the 
5 information that remains to form new tracks. A better approach 
is to integrate both processes where initiating tracks compete 
for the new frame of reports on an equal basis, i.e., at the 
same time. 

In the integrated approach, the data association for track 
10 initiation and track maintenance processing are combined and 
conducted simultaneously. One assignment array is created that 
includes the track scores for potential tracks for both track 
initiation and track maintenance. The first dimension of the 
assignment problem includes only all the tracks either created 
15 or updated in the processing of the prior frames of reports. 
Each of the remaining dimensions accommodates one frame of 
reports . 

After the assignment algorithm finds a solution, each of 
the previously established tracks are updated with the report 
20 assigned to it from the second dimension of the assignment 
array. The remaining reports in the second dimension that were 
in the assignment solution are firmly established as track 
initiators. The unassigned reports in the second dimension of 
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the assignment array are discarded as false signals. The 
processing of the next frame of reports repeats this process 
using these updated tracks and the newly identified initiators 
in the first dimension of the cost array for processing the new 
5 frame of reports. Details of this approach are discussed 
hereinafter . 

The integrated approach just discussed, uses the same 
number of frames of reports for both track initiation and track 
maintenance. The goal in using the multi-dimensional 

10 assignment algorithm is to provide improved performance while 
minimizing the amount of processing required. Typically, track 
initiation will benefit from more frames of reports than will 
track maintenance. ThuSj_ a second approach that integrates 
track initiation and track maintenance is discussed 

15 hereinafter, wherein the number of frames of reports is not the 
same for these two functions. This is a novel approach that is 
introduced in this paper. 

[ ] 

III. 2. 3. Multiple Sensor Processing 

20 There are many advantages and many ways [to] of combining 

data from multiple sensors. There are also many ways of 
categorizing the different algorithmic architectures for 
processing data from multiple sensors. One approach outlines 
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four generic algorithmic architectures (O.E. Drummond. Multiple 
sensor tracking with multiple frame, probabilistic data 
association. In Signal and Data Processing of Small Targets, 
SP1E Proceedings, volume 2561, pages 322-336, 1995) . Two of 
5 these generic architectures are especially pertinent to this 
paper and are summarized briefly. 

In the Centralized Fusion algorithmic architecture, 
reports are combined from the various sensors to form global 
tracks. This algorithmic architecture is also called Central 

10 Level Tracking, Centralized Algorithmic Architecture, or simply 
the Type IV algorithmic architecture. In track maintenance^ 
for example, if a single frame data association is used, then 
for data association a frame of reports from one sensor is 
processed with the latest global tracks; then the global tracks 

15 are updated by the filter function. After the processing of 
this frame of reports is completed, a frame of the reports from 
another (or the same) sensor is processed with these updated 
global tracks. This process is continued as new frames of data 
become available to the system as a whole. 

20 In Centralized Fusion, using the multi-dimensional 

assignment algorithm for the data association with multiple 
sensors is similar to the processing of data from a single 
sensor. Instead of using multiple frames of reports from a 
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single sensor, the multiple frames of reports come from 
multiple sensors. The frames of reports from all the sensors 
are ordered based on the nominal time each frame of reports is 
obtained and independent of which sensor provided the reports. 
5 In this way the approaches discussed in this paper can be 
extended to process reports from multiple sensors using the 
Centralized Fusion algorithmic architecture. 

The second pertinent algorithmic architecture is Track 
Fusion. This approach is also called the Hierarchical or 

10 Federated algorithmic architecture, sensor level tracking or 
simply the Type II algorithmic architecture. In track 
maintenance^ for example, a processor for each sensor computes 
single sensor tracks; these tracks are then forwarded to the 
global tracker to compute global tracks based on data from all 

15 sensors. 

After the first time a sensor processor forwards tracks to 
the global tracker, then subsequent tracks for the same targets 
are cross-correlated with the existing global tracks. This 
track-to-track cross-correlation is due to the common history 
20 of the current sensor tracks and the tracks from the same 
sensor that were forwarded earlier to the global tracker. The 
processing must take this cross-correlation into account and 
there are a number of ways of compensating for this-cross- 
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correlations. One method for dealing with this cross- 
correlation is to decorrelate the sensor tracks that are sent 
to the global tracker. There are a variety of ways to achieve 
this decorrelation and some are summarized in [Drummond, 1996, 
5 Signal Data Proc] recent paper (O.E. Drummond. Feedback in 
track fusion without process noise > In Signal and Data 
Processing of Small Targets [ 19951 , SPIE Proceedings ,[ -2561 : 1 
volume 2561, pages 369-383 , 1995) . 

Once the sensor tracks are decorrelated they can be 

10 processed by the global tracker in almost the same way as 
reports are processed. In the case of track fusion the 
association process is referred to as track (or track-to-track) 
association rather than data (or report-to-track) association. 
If the sensor tracks are decorrelated, the global tracker can 

15 process the tracks from the various sensor processors in much 
the same way that the global tracker of Centralized Fusion 
processes reports. Accordingly the methods described in this 
paper can be readily extended to processing data from multiple 
sensors using either the Centralized Fusion or Track Fusion or 

20 even a hybrid combination of both algorithmic architectures. 

III. 3. Formulation of the Data Association Problem 
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The goal of this section is to briefly outline the 
probabilistic framework for the data association problems 
presented in this work. The technical details are presented 
elsewhere (A.B. Poore. Multidimensional assignment formulation 
5 of data association problems arising from multitarget tracking 
and multisensor data fusion. Computational Optimization and 
Applications, 3:27-57, 1994) . The data association problems 
for multisensor and multitarget tracking considered in this 
work are generally posed (A.B. Poore. Multidimensional 
10 assignment formulation of data association problems arising 
from multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications , 3:21-51, 1994) as 
that of maximizing the posterior probability of the 
surveillance region (given the data) according to 

!5 Maximize {p(T=y | Z N+1 ) \ y 6 r*}, (3.1) 

where Z ml represents N+l data sets, y is a partition of indices 
of the data (and thus induces a partition of the data) , r* is 
the finite collection of all such partitions, r is a discrete 
random element defined on F*, [y° is a reference partition, 

20 ] and P(r=Y|Z N+1 ) is the posterior probability of a partition y 
being true given the data Z N+1 . The term partition is defined 
below; however, this framework is currently sufficiently 
general to cover set packings and coverings. 
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Consider W+l data sets Z{k) ( k=l, . . . , N+l ) , each consisting 
of M k actual reports and a dummy report z£ , and let denote the 
cumulative data set defined by 

5 Z(k) =(z^ =o andZ N+1 = (z(l),... f Z(W+l)} / (3.2) 

respectively. (The dummy report z 0 serves several purposes in 
the representation of missing data, false reports, initiating 

10 tracks, and terminating tracks (A.B. Poore. Multidimensional 
assignment formulation of data association problems arising 
from multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications, 3:27-57, 1994) . ) 
In multisensor data fusion and multitarget tracking the data 

15 sets Z(k) may represent different classes of objects, and each 
data set can arise from different sensors. For track 
initiation the objects are reports that must be partitioned 
into tracks and false alarms. In our formulation of track 
maintenance, which uses a moving window, one data set will be 

20 tracks and remaining data sets will be reports which are 
assigned to existing tracks, as false reports, or to initiating 
tracks . 

We specialize the problem to the case of set partitioning 
(A.B. Poore. Multidimensional assignment formulation of data 
25 association problems arising from multitarget tracking and 

132 



NUMERI .01USC2 

multisensor data fusion. Computational Optimization and 
Applications , 3:27-57, 1994) defined in the following way. 
Define a "track of data" as where each i k [and ] can assume 
rthel zero or nonzero values [ of 0 and , respectively]. A 
5 partition of the data will refer to a collection of tracks of 
data wherein each report occurs exactly once in one of the 
tracks of data and such that all data is used up; the 
occurrence of dummy report is unrestricted. The reference 
partition y° is that in which all reports are declared to be 
10 false. 

Next, under appropriate independence assumptions one can 
show (A.B. Poore. Multidimensional assignment formulation of 
data association problems arising from multitarget tracking and 
multisensor data fusion. Computational Optimization and 
15 Applications, 3:27-57, 1994) 



p ( r- Y |z") n 

p(r= Y °|z N+1 ) (VW 6 Y 

20 i ; , is a likelihood ratio containing probabilities for 
detection, maneuvers, and termination as well as probability 
density functions for report errors, track initiation and 
termination. Define 



25 



*l *N+1 



= -InL. 
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(3.4) 



[1 if]l, if(z. .... ) are assigned to as a track, 

2 1 1 N*l 

[0 otherwise] 0 , otherwise. 



Then, recognizing that 



- In 



P( [ry]r = YlZ Ntl ) 
.p( tv°]r = y°|z n+1 ) 



c. . the problem 



(3.1) can be expressed as the FN+l-dimensional 1 (N+D- 
dimensional assignment problem_ 



Minimize 2. • • • 2, ^.i 



Subject To 



. . _ 1 L \7 i 1 - 1 - 



i =0 i =0 r * N+1 



• * * 2^ zLr • • • z i v ..i H+1 ~ 
ix>r i k = l, ...,M k and £ = 2,...,N, 

2^ * • - 2^ z i r ..i„ ^w+i -1 '-'^!' 



(3.5) 



10 where c ao is arbitrarily defined to be zero. Here, each group 
of sums in the constraints represents the fact that each non- 
dummy report occurs exactly once in a "track of data." One can 
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modify this formulation to include multi-assignments of one, 
some, or all the actual reports. The assignment problem (3.5) 
is changed accordingly. For example, if z ijt is to be assigned 
no more than, exactly, or no less than times, then the " = 1" 
5 in the constraint (3.5) is changed to "<, =, > n i k '" 
respectively. In making these changes, one must pay careful 
attention to the independence assumptions, which need not be 
valid in many applications. Expressions for the likelihood 
ratios can be found in the work of [Poore, 1994, Computational 

10 Optimization & Appl., 3:27-57,1 (A.B. Poore. Multidimensional 
assignments and multitarqet tracking, partitioning data sets, 
I.J. Cox, P. Hansen, and B. Julesz, editors. DIMACS Series in 
Discrete Mathematics and Theoretical Computer Science, American 
Mathematical Society, Providence, R.I., 19:169-198, 1995) and 

15 the references therein. 

III. 4. Track Initiation and Maintenance 
In this section we explain two multiframe assignment 
formulations to the track initiation and maintenance problem. 
20 The continued use of all prior information is computationally 
intensive for tracking, so that a window sliding over the 
frames of reports is used as the framework for track 
maintenance and track initiation within the window. The 
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objectives are to describe an improved version of a simple 
method and then to put this into a more general framework in 
which track initiation and maintenance have different length 
moving windows. 

5 

III. 4.1. The First Approach to Track 
Maintenance and Initiation 

The first method as explained in this section is an 

improved version of our first track maintenance scheme (A. B . 

10 Poore. Multidimensional assignment formulation of data 



association 


problems arisincr from multitaraet tracking 


and 


multisensor 


data fusion. Computational Optimization 


and 



Applications , 3:27-57, 1994) and uses the same window length 
for track initiation and maintenance after the initialization 

15 step. The process is to start with a window of length N+l 
anchored at frame one. In the first step [one] there is only 
one track initiation in that we assume no prior existing 
tracks. In the second and all subsequent frames, there is a 
window of length N anchored at frame k plus a collection of 

20 tracks up to frame k. This window is denoted by {k[; ]j_ 
£+1,..., k+N) . The following explanation of the steps is much 
like mathematical induction in that we explain the first step 
and then step k to step k+1. 

Track Maintenance and Initiation: Step 1. Let 
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{i 1 il 2 ),i 2 [l 2 )^i a il 2 ),i^ 1 U 2 )}^ 



( [4]3. [1]6) 



jbe an enumeration of all those zero-one variables in the 
solution of the assignment problem (3.5) (i.e., z. . =1 ) 



excluding all the false reports in the solution (i.e., all 
those zero-one variables with exactly one nonzero index) and 
zero-one variables in the solution for which (i lf i 2 ) = (0,0). 
(The latter can correspond to tracks that initiate on frames 
three and higher.) These denote our initial tracks. 

Consider only the first two index sets in this enumeration 
([4]3..[1]6) and add the zero index 1 2 =0 with the corresponding 
values of i x and i 2 being zero. Thus, the enumeration is now 



used for this pairing. Suppose now that the next data set f 
i.e., the (N-/-2) th set, is added to the problem. 

To explain the costs for the new problem, one starts with 
the hypothesis that a partition y€T* being true is now 
conditioned on the truth of the pairings on the first two 
frames being correct. Corresponding to the sequence 






the likelihood function is then given by 



(3.7) 




n 



where 



L 
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[4.2a] 

5 Next, define the cost and the corresponding zero-one variable 
C. . . = -lnL 7 . . , 



z i i -i = < 

J -2 i 3 x N+2 



1, if jzi 3 ,..., z^J is assigned to Track T 2 {1 2 ) , 
0, otherwise. 



10 



( [4] 3. [2b] 8) 

respectively. Then the track maintenance problem Maximize 
{ rLylyer*] L v \ yer* \ can be formulated as the following 
multidimensional assignment 



15 



Minimize 



Subject to. 



' " . " . _ J-o-^-a J-mj-O J-o-*"! J- 



1 2 = 0 i 3 =0 i N , 2 =0 



2 3 - t W+2 - L 2-"-3 - t Mt2 



E E 2 ^ =1, 2 2 = 1, [I 2 ,]...,L 2 , 

i 3 =0 i N+2 =0 2 3 »* 2 

E E ••• E z i 3 i,..i m2 =l > i 3 =1 — M 3' 

1 2 =0 i 4 =0 i Nt2 =0 3 3 » +2 

1 2 =0 i 3 =o v^o i p+1 =0 i N+2 =0 2 3 N+2 

for i =1,-...,M and p=4,..., N+2-1, 



(3.9) 



i,=0 



z, , . 6(0,1} for all l 2 ,i 3 ,-,i 
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Track Maintenance and Initiation: Step k. At the beginning 
of the k th step, we solve the following (N+l) -dimensional 
assignment problem. 



Minimize [:]_ 22 22 • • ■ 12 c i ± ± z i ± (3.10) 

Subject to: 



Jt^Jt+l ^Jt+N ^ 



Jf-Jt+I -'-Jt+N 



X v X v • • • X ^ X v • • • X v 



[ Iz,. . = 1[,] 



1 -0 i =() i =0 i =0 i =0 * *+l • ' ' k+N 

for i =1,...,M and p=/c+2, N+Jr-1, 

{4.4) ] 

x k" u 1 k+l" u 1 k+l+N-2~ u 



where for the first step 1 1 and L x are replaced by i x and M lf 
respectively. The first index l k in the subscripts corresponds 
to the seguence of tracks {z^ (l k ) }^* =0 / where 
T k (l k )[={zl^(l k ) r ^z^[l k )}]=^ is a track of 

data from the solution of the problem prior to the formulation 
of ( [4] 3. [4] 10) [ or i] J _If k=lj_ then it is just the first data 
set in frame one. 
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Basic Assumption. Suppose problem ([4] 3. [4] 10) has been 
solved and let the solution, i.e., those zero-one variables 
equal to one, be enumerated by 

(kl^il'^i^Ml-^tV!))}^^ ([4]3.[5]11) 

5 with the following exceptions. 

a. All zero[ ]^one variables for which [l k , i Jt+1 ) = (0, 0) are 
excluded. Thus, tracks that initiate on frames after the 
(k+l) th are not included in the list. 

b. All zero-one variables whose subscripts have the form i k =0 
10 and exactly one nonzero index in the remaining indices 

{i k+1 ,..., i k+N ) are excluded. These correspond to false 
reports on frames rp=k+11 p = Jc+1 , ... , rk+NI Jc+N . 

c. All variables [(]z. . . [)] for which 

(i Jk+1 /.../i Jk+w ) = (0 f 0 f ... / 0) and i k (l k )=0 [in] are excluded from 
15 ( [4] 1. [5]ii) . In other words, the reports on the last AH-1 

frames in the jr. (1. ) , z^ +N , . . . , z\\ are all dummy. Any 
solution with this feature corresponds to a terminated 
track. 

Given the enumeration ( [ 4 ] 3 . [5]UJ , one now fixes the 

20 assignments on the first two index sets in the list 
([4] 3. [5] 11) . The zero index I Jt+1 = 0 is added to the enumeration 
to specify [l k (0) , i Jt+1 (0) )= (0, 0) that is used to represent 
false reports and tracks that initiate on frame k+2 or later, 
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so that the enumeration ( [ 4 ] 3.. [ 5 ] 11. ) is now 

[ (4.6) ] 

Then, for 1. =1, L. . , the lf h such track is denoted by 

AC"*" J. /C + J. Jt+1 

5 T k+1 (l k+1 )={T k (l k (l k ^)) f z^(l k+1 )} and the (N+l)-tuple 

{ T k ^h^> z i2'-"> zk C^\ wil1 denote a track W-W P lus a set 
of reports { z ± k ^ 2 r ' ' r z ^i*n\ f actual or dummy, that are feasible 
with the track T k+1 (l k+l ) • The (W+l) -tuple {^(0), ' Z CLT} 

will denote a track that initiates in the sliding window, i.e., 

10 on subsequent frames. A false report in the sliding window is 
one with all but one non-zero index i p for some p=k+2 , ... ,k+l+N 
in the (N+1) -tuple {^(O), z}£, z£™} . 

The corresponding hypothesis about a partition y 6 T* 
being true is now conditioned on the truth of the L k+1 tracks 

15 existing at the beginning of the N-frame window. (Thus the 
assignments prior to this sliding window are fixed. ) The 
likelihood function is given by 



L = II L i i where 



{ T k+1 U k+1 ) ' Z i k+2 > - ' Z ijt +1+W } 6 Y 



L i i ...i = L t a ) L z^ 2 ... z? +1+N ' 
^i 1 ^ ^ti+N Vi^jc-i' z i k+2 r ' z h+i+N ([4] 3 [7a] 12) 

L W0> =1 - ~~ 
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Next, define the cost and the zero-one variable by 
[said at least one inert gas has an amount, on a molar basis, 
that is greater than the amount of said oxygen in said 
pressurized medium; (4.7b)] 



1 k+l 1 k+2" 1 k+l+N ■ L k*l 1 k*Z 1 k+l+N 1 k+l iJ 'k+l l z l k+2 



0, otherwise, 



i *+l-*Jfc+2"*Jfc+l+W 



(3.13) 



respectively, so that the track extension problem, which was 
originally formulated as Maximize {L y \yer*\, can be expressed as 
exactly the same multi-dimensional assignment in (3.4) [but 
]with k replaced by k+1. Thus, we do not repeat it here. 

10 [ ] 

III. 4.2. The Second Approach to Track 

Maintenance and Initiation 

In the first approach the window lengths were the same for 
both track maintenance and initiation. This can be inefficient 

15 in that one can usually use a shorter window for track 
maintenance than for track initiation. This section addresses 
such a formulation. 
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The General Step k. To formulate a problem for track 
initiation and maintenance we consider a moving window centered 
as frame k of length I+J+l denoted by [k-I r k, k+J] . In this 
representation, the window length for track maintenance is J 
5 and that for initiation, I+J+l. The objective will be to 
explain the situation at this center and then the move to the 
same length window at center k+1 denoted by r rk+1-11 k+1- 

I k+1 rk+l+J1 Jr+1+J 1 , i.e., by moving the r frame 1 window to 

the right one frame at time[ to the right]. The explanation 
10 from the first step follows hereinafter. 

The notation for a track of data is 



T k U*> ={ <U t > ' ~' ' ■■ ([4.8] 3_14 



where the index l k is used for an enumeration of those reports 
paired together. We also use the notation T Pf k^k^ to denote 
15 the sequence of reports belonging to track T k^k^ but 
restricted to frames prior to and including p. Thus, 
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([4]3. [9]15) 

Given this notation for the tracks and partition of the 



data in the frames { k- l f ... f k f ... , k+J} , L T Ptk a k ) will denote the 
accumulated likelihood ratio up to and including frame p(p<k) 
5 for a track that is declared as existing on frame k as a 
solution of the assignment problem. In this notation, the 
likelihood for T Pf k^k^ and that of the association of 

| z i^^"v z i fc *"| w ith track r *(-^) is given by 



L t (i ) =L r a ) L i a j-i a > for an Y P * 
1Q W r^di i vu,) **U*> foranvn< k-l ( [4] ^ 1 [0] & 

h T (1 )i - .V "^r 1 1 ) i 11 >~i 11 )i » i r ° r an Y P s K L r 

respectively. 

The cost for the assignment of j z± k + *r */ j to track 
T k {l k ) and the corresponding zero [ ] -one variable are given by 



c i i ...i = "In [L r ,L. . ] 
25 ^ic+i •'k+j r** 1 *' a Jk+i" ,2 jk+j and 

r J 1, if {z* +1 ,.»/ z* +w I is assigned to track T. (2. ) , 

Z I i i = 1 I J 

* ** J 0, otherwise, 

1 ([4]3.1[1]7) 
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respectively. Likewise, for costs associated with the false 
reports on the frames k-I to k and as associated with 
\ zf +1 ,..., zf +J \ and the corresponding zero-one variables are given 
by: 

ii -i = -ln[L, , . , ] ([4] 3.1 [2] 8) 



Ci k-i"^ 



F 

1 k-l"' 1 k 1 k+l" 1 k+J 



1, if {z ,...,z ,...,z } 

are assigned as a track, 
t 0, otherwise. 



For the window of frames in the range [ [] $k- 1, k [] ] the 
10 easiest explanation of the partition of the reports is based on 
the definitions 

[ (4.13) ] 



T 



i I zf belongs to one of the tracks 
p p . 

listed on frame k, i.e., to one [/ (3*19) 
of the r. (1. ) for some I =1, . . . , I. 



F Pfl =(<i,...,jy\r PiJk )u<o>, 



so that all of those reports not used in the tracks T k {l k ) on 
15 frame p are put into the set of available reports F k for the 
formation of tracks over the entire window. Thus, we formulate 
the assignment problem as 
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Minimize t E S £ c WVw-i M z Lvw-i w 
Subject To £ E £ L<. rVM ,.i w =l O.20) 

{p=W,p*g) i €F ptk r=k*l i r =0 

for i g ef p(t \{0} and q=k-I,...,k, 

£r 52 ^ S Z 2 it . r ..i Jt i^ 1 ...i Jt+J = 1 

for i q =l, ...,M q and q= k+1, k+J, 
z e {0,1} for all i ,...,i ,i ,...,i , 

1 k-r 1 k 1 k*l" 1 k+J * 1 K K * L K ° 

£ 4.1-^ = 1 ' 2 * = 1 * - * v 

r~k+l i r =0 

/I ^1 z l u i M „i i CT =1 '-' 



l k =l {r=k+l, r*q) i r =0 

for q=k+l, k+ J, 
z n i E {°^ 1 } for a11 l-k* 1 ' ik+i'-' ik+j- 

1 k 1 k+l" 1 k+J K K x K ° 



[Minimize c Ic . I „ iIC(ik+1 „. iIC+J z k . I „ ikik+1 ... ik+J Subject to: for 
irl,-,M k and k=2,...,N, zii"l N+ i ^{0,1} 

for all ± lt • ~, i N+1 ] 
5 Basic Assumptions. Suppose problem ( T4 . 141 1.20 ) has been 

solved and let the solution, i.e.[,] those zero-one variables 
equal to one, be enumerated by 

, k ^: L ([4] 3. [15121) 

\z F . r +1 

with the following exceptions. 
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a. All zero [ ]^one variables in the second list for which 
(i k _ x ,...,i k ,i k+1 ) = [0,...,0,0) are excluded. Thus, tracks that 
initiate on frames after the (Jc+1) th are not included in 
the list. 

5 b. All false reports are excluded, i.e., all zero-one 

variables in the second list whose subscripts have exactly 

one nonzero index, 
c. All variables z/, , for which (i*. +1 , i t+M ) = (0, 0, 0) and 

Z[p) DT k {l k ) =0 for p=k-K,...,k where \KzO MC is user 

10 specified. Thus the track T k [l k ) is terminated if it is 

not observed over K+J+l frames. 
Given the enumeration ( [4] 3. [15] 20) , one now fixes the 
assignments on the all index sets up to and including the 
(k+l) th index sets. 



T 11 ) = 4 



f ° r 1 k+l = 1 ' ■ ' "\ tl > 

for 1 *+l =£ k.l + 1 '-" L *+l- 



(3.22) 



15 [ (4.16)] 
Thus one can now formulate the assignment problem for the 
next problem exactly as in ( [4 . 141 3.20 ) but with k replaced by 
k+1. Thus, we do not repeat it here. 
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The Initial Step. Here is one explanation for the initial 
step. First, assume that N=I+J. In this case, we start the 
track initiation with a solution of (3.5) . Let 
{ i, ( I 2 ) , i 2 ( 1 2 ) , - , i N ( i 2 ) , i ml ( 1 2 ) }^ =Q 
5 be an enumeration of the solution set of (3.5), i.e., those 
zero-one variables z i l( i 2 ,i 2 (i 2 ,..-i N+lU2 ) =1 , including [z 00 ... 0 =l] z OQ ... 0 =\ 
corresponding to 12=0, but excluding all those zero-one 
variables that are assigned to one and correspond to false 
reports (i.e., there is exactly one nonzero index in the 
10 subscript of z ± . . ) a ll those zero-one variables that are 

r -4-4 2 w+i ' 

assigned to one and correspond to tracks that initiate on 
frames higher than J+2. Then we fix the data association 
decisions corresponding to the reports in our list of tracks 
prior to and including frame Jc+l=X+2. This defines the k for 
15 problem ([4] .3. [4] 5) and one can then continue the development 
by adding a frame to the window as in the general case . 

If [ I+J>N] J+J" > N , then one possibility is to start the 
process with frames, and assuming J<N, proceed as before 

replacing I by N-J for the moment, and continue to add frames 
20 without lopping off the first frame in the window until reaches 
a window of length I+J+l. Then we proceed as in the previous 
paragraph . 

If [ I+J<N] J+J" < N , then one can solve the track initiation 
problem (3.5), formulate the problem with the center of the 
25 window at [Jc+l=Z\r+l - Jl k+1 = N+l-J , enumerate the solutions as 
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above, and lop off the first N-J-I frames. Then, we proceed 
just as in the case I+J=N. 

III. 5. CONCLUDING COMMENTS 

A primary objective in this work has been to demonstrate 
5 how multidimensional assignment problems arise in the tracking 
environment. The problem of track initiation and maintenance 
has been formulated within the framework of a moving window 
over the frames of data. The solution of these NP-hard, noisy, 
large scale, and sparse problems to the noise level in the 

10 problem is fundamental to superior track estimation and 
identification. Thus, one must utilize the special structure 
in the problems as well take advantage of special information 
that is available. Since these moving windows are overlapping, 
there are some algorithm efficiencies that can be identified 

15 and that take advantages of the overlap in the windows from one 
frame of reports to the next. Here is an example of the use of 
a primal solution of one problem to warm start the solution of 
the next problem in the sequence. 

Suppose we have solved problem ([4] 3. [4] 10) and have 

20 enumerated all those zero-one variables in the solution of 
([4] 3. [4] 10) as in ( [4] 3 . [5] 11) . Add the zero index l k+1 =0 , so 
that the enumeration is 

With this enumeration one can define the cost by 
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c i i =c i(i nil (l )i i i (l )i ([512.24) 



and the two [ ] ^dimensional assignment problem 



^ 0 = Minimize J* !C c i i z i i = V 2 (z 2 ) 

2 Jt+i = u 1 Jk+i+w =u 



Subject to X Z v w+N =1 ^ 2 * + i = 1 '-' L ;c + i' 



2 z i +1 i Jk+1+lf 1 ' 2 jk+i+w 

Z 2 w i W ^ |0 ' 1} 311 2 **1' 



[(5.3)] 



(3.25) 



Let w be an optimal or feasible solution to this two- 
dimensional assignment problem and define 



10 



z . 



" 1 k*H 1 k*l 



and ^ 7 , =1 for some 2. + =1, L. +1 , p , 

* +1 * +1 ([512- [4] 26) 

or if (I Jt+1 .i Jt+1+w ) = (0,0) 
0, otherwise. 



This need not satisfy the constraints in that there are usually- 
many objects left unassigned. Thus, one can complete the 
15 assignment by using the zero-one variables in ( [4] 2- [4] 10) with 
k replaced by Jc+1 with exactly one nonzero index corresponding 
to any unassigned object or data report. 
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For the dual solutions, the multipliers arising from the 
solution of the two [ ] ^dimensional assignment problem 
( rs . 31 3 . 25 ) corresponding to the second variable, i.e., 
(u k+1+N P +1+w . These are good initial values [for] to use in 
5 a relaxation scheme r Til . 121 1 (A.B. Poore an d N.Riiavec. 
Partitioning multiple data sets: multidimensional assignments 
and Lagrangian relaxation, in quadratic assignment an d related 
problems; P.M. Pardalos and H. Wolkowicz, editors , DIMACS 
series in Discrete Mathematics and Theoretical Computer 

10 Science. 16:25-37, 1994: A.J. Robertson III. A class of 
lagrangian relaxation algorithms for the multidimensional 
assignment problem. Ph.D. Thesis, Colorado State University, 
Ft. Collins, CO, 1995) . Finally, note that one can also 
develop a warm start for problem ( [4] 3. [14] 20.) in a similar 

15 fashion. 



IV. REVIEW OF NEW RELAXATION SCHEMES 

IV. 1. Introduction 

Multidimensional assignment problems govern the central 

20 problem of data association in multisensor and multitarget 
tracking, i.e., the problem of partitioning observations from 
multiple scans of the surveillance region and from either 
single or multiple sensors into tracks and false alarms. This 
fundamental problem can be stated as (A.B. Poore . 

25 Multidimensional assignments and multitarget tracking. 
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partitioning data sets, I.J. Cox, P. Hansen, and B. Julesz, 
editors, DIMACS Series in Discrete Mathematics and Theoretical 
Computer Science, American Mathematical Society, Providence, 
R.I. , 19:169-198, 1995; A. B . Poore . Multidimensional assignment 
5 formulation of data association problems arising from 
multitarget tracking and multisensor data fusion. Computational 
Optimization and Applications r 3:27-57, 1994) 



Minimize 




N 



(4.1) 



Subject to: 




[(1.1) 




for i p =l,...,M p and p=2,..., 



N-l, 




for i k =l,...,M k and k=2,... ,N, 



152 



NUMERI .01 USC 2 

zii-l^x e{0,l} for all t xt - t L^ ] 

where c 0 ... 0 is arbitrarily defined to be zero and is included for 
notational convenience. One can modify this formulation to 
5 include fmultiassignmentsl mult i -assignments of one, some, or 
all of the actual reports. The zero index is used in 
representing missing data, false alarms, initiating and 
terminating tracks as described in (A.B. Poore. 
Multidimensional assignments and multitarget tracking, 

10 partitioning data sets, I.J. Cox, P. Hansen, and B. Julesz, 
editors, DIMACS Series in Discrete Mathematics and Theoretical 
Computer Science, American Mathematical Society, Providence, 
R.I., 19:169-198, 1995; A.B. Poore. Multidimensional assignment 
formulation of data association problems arising from 

15 multitarget tracking and multisensor data fusion. Computational 
Optimization and Applications, 3:27-57, 1994) . In these 
problems, we assume that all zero-one variables z- ■ with 

J -1"' J -N 

precisely one nonzero index are free to be assigned and that 
the corresponding cost coefficients are well-defined. (This is 

20 a valid assumption in the tracking environment (A.B. Poore. 
Multidimensional assignments and multitarget tracking, 
partitioning data sets, I.J. Cox, P. Hansen, and B. Julesz, 
editors, DIMACS Series in Discrete Mathematics and Theoretical 
Computer Science, American Mathematical Society, Providence , 

25 R.I., 19:169-198, 1995; A.B. Poore. Multidimensional assignment 
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formulation of data association problems arising from 
multitarget tracking and multisensor data fusion. Computational 
Optimization and Applications, 3:27-57, 1994) .) Although not 
required, these cost coefficients with exactly one nonzero 
5 index can be translated to zero by cost shiftin g (A. B . Poore 
and N. Rijavec. A lagrangian relaxation algorithm for 
multidimensional assignment problems arising form multitarget 
tracking. SI AM Journal of Optimization, 3, No. 3:544-563, 1993) 
without changing the optimal assignment. Finally, this 
10 formulation is of sufficient generality to include the 
symmetric problem and the asymmetric inequality problem [.] 
(A.J. Robertson III. A class of lagrangian algorithms for the 
multidimensional assignment problem. Ph.D. Thesis, Colorado 
State University, Ft. Collins, CO, 1995) . 

15 

The data association problems in tracking that are 
formulated as Equation ( Til 4 . 1) have several characteristics. 
They are normally sparse, the cost coefficients are noisy and 
the problem is NP-har d (M.R. Galey and D.S. Johnson. Computers 

20 and Intractability. W.H. Freeman and Company, San Francisco, 
CA, 1979) , but it must be solved in real-time. The only known 
methods for solving this NP-hard problem optimally are 
enumerative in nature, with branch- and-bound being the most 
efficient; however, such methods are much too slow for real- 

25 time applications. Thus one must resort to suboptimal 
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approaches. Ideally, any such algorithm should solve the 
problem to within the noise level, assuming, of course, that 
one can measure this noise level in the physical problem and 
that the mathematical method provides a way to decide if the 
5 criterion has been met. 

There are many algorithms that can be used to construct 
suboptimal solutions to NP-hard combinatorial optimization 
problems. These include greedy (and its many variants), 

10 relaxation, simulated annealing, tabu search, genetic 
algorithms, and neural [netork algorithms! network algorithms 

(C.H. Papadimitriou and K.Steiglitz. Cowbi na t ori a 1 

Optimization: Algorithm and Complexity. Prentice-Hall, Inc., 
Engiewood Cliffs, NJ, 1982; J. Pearl. Heuristics: Intelligent 

15 Search Strategies for Computer Problem Solving. Addison-Wesley , 
Reading, MA, 1984; C.R. Reeves ed. Modern Heuristic Techniques 
for Combinatorial Problems. Halstead Press, Wiley, NY, 1993) . 
For the three^ dimensional assignment problem, Pierskalla (W. 
Pierskalla. The tri-substitution method for three-dimensional 

20 assignment problem. Journal du CORS, 5:71-81, 1967) developed 
the tri-substitution method, which is a variant of the simplex 
method. Frieze and Yadegar (A.M. Frieze and J. Yadegar. An 
algorithm for solving 3 -dimensional assignment problems with 
application to scheduling a teaching practice. Journal of the 

25 Operational Research Society, 39:989-955, 1981) introduced a 
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method based on Lagrangian relaxation in which a feasible 
solution is recovered using information provided by the relaxed 
solution. Our choice of approaches is strongly influenced by 
the need to balance real-time performance and solution quality. 
5 Lagrangian relaxation based methods have been used successfully 
in prior tracking applications T . 1 (K.R. Pattipati, S. Deb, and 
Y. Bar-Shalom. A s-dimensional assignment algorithm for track 
initiation. In Proceedings of the IEEE Systems Conference, 
Kobe, Japan, pages 127-130, 1992; Y. Bar-Shalom, S. Deb, K.R. 

10 Pattipati, and H. Tsanakis. A new algorithm for the generalized 
multidimensional assignment problem. In Proceedings of the IEEE 
International Conference on Systems, Math, and Cybernetics, 
Chicago, pages 132-136, 1992; A.B. Poore and N. Riiavec. A 
numerical study of some data association problems arising in 

15 multitarget tracking. Large Scale Optimization: State of the 
Art, W.W. Hager, D.W. Hearn and P.M. Pardalos, editors. Kluwer 
Academic Publishers B.V. . Boston, pages 339-361, 1994; A.B. 
Poore and N. Riiavec. Partitioning multiple data sets: 
multidimensional assignments and lagrangian relaxation. In P.M. 

20 Pardalos and H. Wolkowicz, editors. Quadratic assignment and 
related problems: DIMACS Series in Discrete Mathematics and 
Theoretical Computer Science, volume 16, pages 25-37, 1994; 
A.B. Poore and N. Riiavec. A lagrangian relaxation algorithm 
for multidimensional assignment problems from multitarget 

25 tracking. SIAM Journal of Optimization. 3 . No. 3:544-563, 1993). 
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An advantage of these methods is that they provide both an 
upper and lower bound on the optimal solution, which can then 
be used to measure the solution quality. These works extend 
the method of Frieze and Yadegar (A.M. Frieze and J.Yadegar. An 
5 algorithm for solving 3 -dimensional assignment problems with 
application to scheduling a teaching practice. Journal of the 
Operational Research Society. 32:989-995, 1981) to the 
multidimensional case . 

10 IV. 2. Probabilistic Framework for Data Association. \ (ABP) 1 

The goal of this section is to explain the formulation of 
the data association problems that governs large classes of 
data association problems in centralized or hybrid centralized- 

15 sensor level multisensor/multitarget tracking. The 
presentation is brief; technical details are presented for both 
track initiation and maintenance in the work of rPoorel (A.B. 
Poore. Multidimensional assignments and multitarget tracking, 
partitioning data sets, I.J. Cox, P. Hansen, and B. Julesz, 

20 editors, DIMACS Series in Discrete Mathematics and Theoretical 
Computer Science, American Mathematical Society, Providence, 
R.I., 19:169-198, 1995; A.B. Poore . Multidimensional assignment 
formulation of data association problems arising from 
multitarget tracking and multisensor data fusion. Computational 

25 Optimization and Applications, 3:27-57, 1994) . The formulation 
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is of sufficient generality to cover the MHT work of Reid, 
Blackman and Stein (S.S. Blackman. Multiple Target Tracking 
with Radar Applications . Artech House, Norwood, MA, 1986) and 
modifications by Kurie n (T.Kurien. Issues in the designing of 
5 practical multitarget tracking algorithms. In Mul ti target- 
Multisensor Tracking: Advanced Applications by Y. Bar-Shalom. 
Artech House . MA, 1990) to include maneuvering targets. As 
suggested by Blackma n (S.S. Blackman. Multiple Target Tracking 
with Radar Applications, Artech House, Norwood, MA, 1986) , this 

10 formulation can also be modified to include target features 
(e.g., size and type) into the scoring function. The recent 
work [of Poore and Drummondl (A.B. Poore and O.E. Drummond. 
Track initiation and maintenance using multidimensional 
assignment problems. In D.W. Hearn, W.W. Hager, and P.M. 

15 Pardalos, editors. Network Optimization, volume 450, pages 407- 
422, Boston, 1996. Kluwer Academic Publishers B.V.) 
significantly extends the work of this section to new 
approaches for multiple sensor centralized tracking. Future 
work will involve extensions to track-to-track correlation. 

20 

The data association problems for multisensor and 
multitarget tracking considered in this work are generally 
posed as that of maximizing the posterior probability of the 
surveillance region (given the data) according to 
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(4^.2 [.1] ) 



where Z N represents N data sets, y is a partition of indices of 
the data (and thus induces a partition of the data) , r* is the 
finite collection of all such partitions, r is a discrete 
5 random element defined on r*, [y° is a reference partition,] and 
P(rr=y|Z N l r = y \ Z N ) is the posterior probability of a partition 
Y being true given the data Z N . The term partition is defined 
below; however, this framework is currently sufficiently 
general to cover set packings and coverings (A.B. Poore. 

10 Multidimensional assignment formulation of data association 
problems arising from multitarget tracking and multisensor data 
fusion. Computational Optimization and Applications, 3:21-51, 
1994; A.B. Poore. Multidimensional assignments and multitarget 
tracking, partitioning data sets, I.J, cox, P. Hansen, and B. 

15 Julesz, editors, DIMACS Series in Discrete Mathematics and 
Theoretical Computer Science, American Mathematical Society, 
Providence, R.I., 19:169-198, 1995). 



Consider N data sets Z{k) {k=l, ...,N) each of M k reports 



20 




and let Z N denote the cumulative data set defined by 
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( [2]4. [2] 3) 

respectively. In multisensor data fusion and multitarget 
tracking the data sets Z{k) may represent different classes of 
objects, and each data set can arise from different sensors. 
5 For track initiation the objects are measurements that must be 
partitioned into tracks and false alarms. In our formulation 
of track maintenance (A.B. Poore. Multidimensional assignment 
formulation of data association problems arising from 
multitarget tracking and multisensor data fusion. Computational 

10 Optimization and Applications, 3:27-57, 1994; A.B. Poore , 
Multidimensional assignments and multitarget tracking, 
partitioning data sets, I.J, cox, P. Hansen, and B. Julesz, 
editors, DIMACS Series in Discrete Mathematics and Theoretical 
Computer Science, American Mathematical Society* Providence, 

15 R.I., 19:169-198, 1995) , which uses a moving window over time, 
one data set will be tracks and remaining data sets will be 
measurements which are assigned to existing tracks, as false 
measurements, or to initiating tracks. In sensor level 
tracking, the objects to be fused are tracks (S.S. Blackman. 

20 Multiple Target Tracking with Radar Applications . Artech House, 
Norwood, MA, 1986) . In centralized fusio n (S.S. Blackman. 
Multiple Target Tracking with Radar Applications . Artech House, 
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and Z w = {z(l),...,Z(N)}, 
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Norwood, MA, 1986) , the objects may all be measurements that 
represent targets or false reports, and the problem is to 
determine which measurements emanate from a common source. 

5 We specialize the problem to the case of set partitioning 



(A.B. Poore. 


Multidimensional assicrnment formulation of data 


association 


problems arisincr from multitaraet trackina and 


multisensor 


data fusion. Computational Optimization and 


Applications 


3:27-57, 1994; A.B. Poore . Multidimensional 


assignments 


and multitaraet trackina, partitionina data sets. 


I.J. cox, P. 


Hansen, and B. Julesz. editors, DIMACS Series in 


Discrete Mathematics and Theoretical Computer Science, American 


Ma thema tical 


Society, Providence, R.I.. 19:169-198, 1995) 


defined in 


the following way. First, for notational 



5 convenience in representing tracks, we add a zero index to each 
of the index sets in Equation ( [2] 4. [2] 3) , a dummy report z k 0 to 
each of the data sets Z(k) in Equation ( [2] 4. [2] 3) , and define 

a "track of data" as ( z v~' z £ ) ' where i k can now assume the 
value [s] of 0[ and respectively] . A partition of the data will 

0 refer to a collection of tracks of data wherein each report 
occurs exactly once in one of the tracks of data and such that 
all data is used up; the occurrence of dummy report is 
unrestricted. The dummy report z 0 * serves several purposes in 
the representation of missing data, false reports, initiating 

5 tracks, and terminating tracks (A.B. Poore. Multidimensional 
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assignment formulation of data association problems arising 
from multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications, 3:27 -57, 1994; 
A.B. Poore. Multidimensional assignments and multitarget 
5 tracking, partitioning data sets, I.J, cox, P. Hansen, and B. 
Julesz, editors, DIMACS Series in Discrete Mathematics and 
Theoretical Computer Science, American Mathematical Society, 
Providence, R.I.. 19:169-198, 1995) . The reference partition 
is that in which all reports are declared to be false. 

10 

Next under appropriate independence assumptions one can 
show (A.B. Poore. Multidimensional assignment formulation of 
data association problems arising from multitarget tracking and 
multisensor data fusion. Computational Optimization and 
15 Applications, 3:21-51, 1994; A.B. Poore. Multidimensional 
assignments and multitarget tracking, partitioning data sets, 
I.J, cox, P. Hansen, and B. Julesz, editors, DIMACS Series in 
Discrete Mathematics and Theoretical Computer Science, American 
Mathematical Society, Providence, R.I.. 19:169-198, 1995) 

20 P(r = Y l ZN) ^ L = , . n v L. . , ( [2]4. [3]4) 

p(r= Y °|z N ) Y (V 2 W > e Y ^ 



where y° is a reference partition, L i r ..i N is a likelihood ratio 
containing probabilities for detection, maneuvers, and 
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termination as well as probability density functions for 
measurement errors, track initiation and termination. Then if 
c. . = -InL. 



-In 



P(r=y|Z N ) 
P(r = Y°|Z N ) . 



( [2]4. [4] 5) 



Using Equation ( [2] 4 . [33 4) and the zero-one variable 
z. . = 1 if (i-,-,i w ) 6 y and 0 otherwise, one can then write the 
10 problem ([2]4.[1]4) as the following ^-dimensional assignment 
problem: [ (2.5) ] 



Minimize \\ c. . z. . 

Subject To £ z i i =1 U^l,-,^) , (4.6) 

i 2 i 3 ..-i„ 

. £ Z Vi„ =1 (i P = 1 '" ' M p axid P = 2,...,N-1) , 

i l i 2*" i N-l 

z. .6 {0,1} for all i , 

where c 0 ... 0 is arbitrarily defined to be zero. Here, each group 
of sums in the constraints represents the fact that each-non- 
15 dummy report occurs exactly once in a "track of data." One can 
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modify this formulation to include rmultiassignmentsl multi- 
assignments of one, some, or all the actual reports. The 
assignment problem Equation ( [2] 4. [5] 6) is changed accordingly. 
For example, if z^ is to be assigned no more than, exactly, 
5 or no less than n± times, then the w = l" in the constraint 
([2]4. [5]6) is changed to <, =, >n£, !l respectively. 

Modifications for group tracking and multi^resolution features 
of the surveillance region will be addressed in future work. 
In making these changes, one must pay careful attention to the 
10 independence assumptions that need not be valid in many 
applications . 

Expressions for the likelihood ratios -^i ,i [,] can be 
found in the work of Poore (A.B. Poore . Multidimensional 
assignment formulation of data association problems arising 

15 from multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications, 3:21-51. 1994; 
A.B. Poore. Multidimensional assignments and multitarget 
tracking, partitioning data sets, I.J, cox, P. Hansen, and B. 
Julesz, editors, DIMACS Series in Discrete Mathematics and 

20 Theoretical Computer Science, American Mathematical Society, 
Providence, R.I., 19:169-198, 1995) . These expressions include 
the developments of Reid, Kurien (T. Kurien. Issues in the 
designing of practical multitarget tracking algorithms. In 
Multitaraet-Multisensor Tracking: Advanced Applications by Y. 

25 Bar-Shalom. Artech House, MA, 1990) , and Stein and Blackman 
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(S.S. Blackman. Multiple Target Tracking with Radar 
Applications. Artech House, Norwood, MA, 1986) . What 1 s more , 
they are easily modified to include target features and to 
account for different sensor types. In tracJc initiation, the 
5 N data sets all represent reports from N sensors, possibly all 
the same. For track maintenance, we use a sliding window of N 
data sets and one data set containing established tracks. (A.B . 
Poore. Multidimensional assignment formulation of data 
association problems arising from multitarget tracking and 

10 multisensor data fusion. Computational Optimization and 
Applications, 3:27-57, 1994; A.B. Poore. Multidimensional 
assignments and multitarget tracking, partitioning data sets, 
I.J, cox, P. Hansen, and B. Julesz, editors, DIMACS Series in 
Discrete Mathematics and Theoretical Computer Science, American 

15 Mathematical Society, Providence, R.I.. 19:169-198, 1995). The 
formulation is the same as above except that the dimension of 
the assignment problem is now AT+1. [@@@] 

[ To explain the costs for the new problem, one starts with 
the hypothesis that a partition y e r* being true is now 
20 conditioned on the truth of the pairings on the first two 
frames being correct. Corresponding to the sequence 

{T 0 (l 0 ) , z. z. }, the likelihood function is then given by 

1 <• 2 3'"' 2 N+2 
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n L, . »• . where 

iT 2 (l 2 ) , z ,...,z. } ey 

L 7 . = L . j . , z . , ...Z. = L {1 .L t ...z. t 

L m = 1. 



Next, define the cost and the corresponding zero-one variable 

(4 . 2b) 



'i i -i = " lni ii -i ' 

J -2 i 3 ^ilM ■ 4 2 i 3 i N+2 



.-J 3 N + 2, 
if \z L ,...,z. \ . 

1 -"-3 i w+2 j 



i 2 i,..i Nt2 " | 1 " I'V-'^i^f is assigned to Track r 2 (J 2 ), 



0 otherwise 



respectively. Then the track maintenance problem Maximize 
5 {Ly|Y Er *} can be formulated as the following multidimensional 
assignment 



Minimize 



5 s f Y c 



1,-0 i 3 =0 i Nt2 =0 



2* t 3"' J, N 2^3' ■ L N+2 
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Subject to: 

A " " i -0 ~ 1 ' 12 1 ' "^2 ' ' • - ' i 2' 

1 3~ U 1 N+2~ U 

2^ 2^ • - • 2^ z i,i v ..i w+? 1 ^ 2 3 1/ * • • * M 3' 

l=0i=0 i =0 33 W+2 

2^ 2^ • • • 2^ 2^ • • • 2-r z i,i,».i„„ 1 



1 2 =0 i 3 = 0 i p+1 = 0 i N+2 = 0 2 3 " +2 



for i k =l,...,M k and k=2,...,N, 



ii=0 



zi r -l N+1 6(0,1} for all i 1# -,i N+1 

5 

Track Maintence and Initiation: Step k. At the beginning 
of the k th step, we solve the following (N+l) -dimensional 
assignment problem. 



Minimize 



/ ' / J • • • / ' ° 7 ... 7 ... i 

k k+1 k+N 
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Subject to: 

2^ — .2^. z i,i,...i ~ 1 ' - i 2 -1 ' i 2 /...,£ 2 / 



V 0 W> 



1-0 i -o ' * 'i -o ^VV." 1, ^s" 1, * ' * ' M3 ' 



V 0 i 3 =0 Vl = ° Vl = ° W 0 



for i p =l,...,Mp and p=k+2 , ... , N, +k+l , (4.4) 




zl k i k+1 -l k+N e{0,l} for all l k/ i k+1/ i k+N 

where for the first step l x and L x are replaced by i x and M lf 

respectively. The first index l k in the subscripts corresponds 

r A 

to the sequence of tracks \T k {l k )!^ = 0 where 

J 1 k \ k 

T k^k^ ~^ z i ( ) / "> z ^ (1^' a track °f data from the solution 

of the problem prior to the formulation of (4.4) or if k=l then 
it is just the first data set in frame one. 
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Suppose problem (4.4) has been solved and let the 
solution, i.e., those zero-one variables equal to one, be 
enumerated by 

{ (l k (l k+1 ) ,i k+1 (l k+1 ) ,...,i k+ Ndk + i) ) }L k+1 l k+1 =l (4.5) 
5 with the following exceptions. 

a. All zero one variables for which (l k , i k+1 ) = (0, 0) are 
excluded. Thus, tracks that initiate on frames after the 
(k+l) th are not included in the list. 

b. All zero-one variables whose subscripts have the form l k =0 
10 and exactly one nonzero index in the remaining indices 

{i k+1 , i k+N } are excluded. These correspond to false 
reports on frames p=k+l, k+N. 

c. All variables z lkik+1 „ ik+N ) for which (i k+1 ,..., i k+N ) = (0, 0, . . . , 0) 
and i k (l k )=0 in (4.5) . In other words, the reports on the 

15 last N+l frames in the )T. (1.), ^ | are all dummy. Any 

solution with this feature corresponds to a terminated 
track. 

Given the enumeration (4.5), one now fixes the assignments on 
the first two index sets in the list (4.5). The zero index 
20 l k+ i=0 is added to the enumeration to specify 
T k (l k ) ={z\ (l k ),...,zl U k )) that is used to represent false 
reports and tracks that initiate on frame k+2 or later, so that 
the enumeration (4.5) is now 
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Then, for 2 Jt+1 = 1, L k+1 , the such track is denoted by 

W^+iMvV^+i))' and the (N+l) -tuple 

/ k+2 Jt+l+AH + , 

\T {1 ) , z i ,...,Zj t will denote a track Z\ +1 (l. J plus a 
v K+i *+± ± ,*+i+n J , ^ 

/ k+Z k+l+Nl 

set of reports \z i f - 9 Zx f, actual or dummy, that are 

1 J_ k+2 X k+1+N ; 

feasible with the track T k+1 (l k+1 ) . The (N+l) -tuple 

f k+2 k+l+Nl 

IT. ^O), za ," f z\ f, will denote a track that initiates m 

1 * 1 x k+2 X k+1+N } 

the sliding window, i.e., on subsequent frames. A false report 
in the sliding window is one with all but one non-zero index i p 

for some p = k + 2 , ... , k + l+N in the (N+l) -tuple 

/ rt k+2 k+l+Nl 
\T k+1 (0), z ± ,-, Zi ). 

1 K 1 x k+2 "'"k+l+N } 

The corresponding hypothesis about a partition y e r* being 
true is now conditioned on the truth of the L k+1 tracks existing 
at the beginning of the N- frame window. (Thus the assignments 
prior to this sliding window are fixed.) The likelihood 
function is given by 



n 

T k (1 )' Z L ' •' Z i e Y 



L = 11 L. . t where 



(4.7a) 

= r. r . o ... t. 



L, . . - L T n v f L k+2/"'/i k+l+N 



Jt+l- L Jt+2"* J -Jt+l+N A Jt+l * A Jt+l ' Zj Z. 

■ L k+2 """k+l+N 
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Next, define the cost and the zero-one variable by 

, k+2 k+l+N 
Z 7 . = f 1 if \z j_ z ± 

x *+i 2 Jt+i+w J k+2 k+l+N 

I 0 otherwise, (4.7b) 

respectively, so that the track extension problem, which was 
5 originally formulated as Maximize {L^yer*}, can be expressed as 
exactly the same multi -dimensional assignment in (3.4) but with 
k replaced by k+1. Thus, we do not repeat it here. 

The Second Approach to Track Maintenance and Initiation 
In the first approach the window lengths were the same for 
10 both track maintenance and initiation. This can be inefficient 
in that one can usually use a shorter window for track 
maintenance than for track initiation. This section addresses 
such a formulation. 

The General Step k. To formulate a problem for track 
15 initiation and maintenance we consider a moving window centered 
as frame k of length I+J+l denoted by [k-I,...,k,...,k+J] . In this 
representation, the window length for track maintenance is J 
and that for initiation, I+J+l. The objective will be to 
explain the situation at this center and then the move to the 
20 same length window at center k+1 denoted by [k+1- 
I,...,k+l,...,k+l+J] , i.e., by moving the frame one frame at time 
to the right. The explanation from the first step follows 
hereinafter. 

The notation for a track of data is 
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W=< z i ,<V'"'*i ^ k )>",zl (!,)> (4-8) 



where the index l k is used for an enumeration of those reports 
paired together. We also use the notation T pk (l k ) to denote 
the sequence of reports belonging to track T k (l k ) but restricted 
5 to frames prior to and including p. Thus, 



l P 



(4.9) 



Given this notation for the tracks and partition of the data in 
the frames {k-I, k, k+j} , LT pk (l k ) will denote the accumulated 
likelihood ratio up to and including frame p(p^k) for a track 
10 that is declared as existing on frame k as a solution of the 

assignment problem. In this notation, the likelihood for 

. . / k+1 k+N\ . , 

T nkdk) and that of the association of \z i ,-,Zi f with track 

T k {l k ) is given by 

L Tii)i ...i =L r { i) L i iD-iii )± foranyp*k-l, 



172 



NUMERI .01 USC 2 

respectively. 

^ / k+1 k+N\ 
The cost for the assignment of \z ± r" ,?>\ \ to track 

1 X k+1 X k+N' 

T k {l k ) and the corresponding zero one variable are given by 



1 i -i = " ln [ L r (i ) L i -i ] 

Jc k+1 k+J 1 Jt^Jt' 1 k^i ± k*j 



and 



= -ln[L 



x * x *+l X *+J 



J k + 1 k+N\ . . ,., % 

11 z H z j is assigned to track T. (1. ) , 

l x k-l x k-l/ 

0 otherwise, 



(4.11) 

respectively. Likewise, for costs associated with the false 

reports on the frames k-I to k and as associated with 

Jt+1 k+J\ ^ ^ JB . ^ 

z 7 * , ...,z 7 - f and the corresponding zero-one variable are given 



by: cFik.j »-i k i k+1 -i k+J = -ln[l, J 

L k-I k k+1 k+J J 



Z Fl k _i"lklk + l"lk+J - 



1 { if z. k i ,...,z. ki ,...,zi k+J \ 

are assigned as a track, 
0 otherwise. 



(4.12) 



For the window of frames in the range [k-I,...,k], the 
easiest explanation of the partition of the reports is based on 
the definitions 



F - < 



• p belongs to one of the tracks listed on frame k, 
x p lz i i.e., to one of the T k {l k ) for some l k =l,...,L k 



= ({!,..., Mp)\T ptk ) u(0} 



(4.13) 
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so that all of those reports not used in the tracks T k {l k ) on 
frame p are put into the set of available reports F pk for the 
formation of tracks over the entire window. Thus, we formulate 
the assignment problem as 

5 



V V* F F 

7 j 7 j 7 v 7 j C 7 7 -1 7' Z 7* 7* 7* -j 

P =k-i i p €F P(k r=k+i i t =o ± k-r-- L k' ■ L k+i~- L k+j ■ L k-r"- L k ± k+v- L k 

^ w f r T 



_ _ . ..... _ _ k+J 

Minimize 



l k =l r=k+l i r =0 ^k-k+V"-k+J ^k^k+V^k+J 



Subject to: 



1, r '"^k^k+l" k+J 

t~= IrA- "I i =(1 

P P/J 

for i-6F_A{0} and q=k-I,...k, 



^k-I,p*q) i p eF p/Jt r=*+l i r =0 



g g,/c 
=Jt-I i p eF pJt {r=*+l,r#gl i r =0 



r WV^* +J 6{0 ' 1} for a11 ik-T'-'ik'^' -'h 

£ Z l,i, .-i, = If i,= lr~f*V 

=Jc+l i r =0 * ^/c+J 9 9 

for q= k+1, k+J, 
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solution, i.e., those zero-one variables equal to one, be 
enumerated by 




(4.15) 



5 with the following exceptions. 

a. All zero one variables in the second list for which 
(i w ,..., i k r i k+1 ) = (0,...,0,0) are excluded. Thus, tracks that 
initiate on frames after the (k+l) th are not included in 
the list. 

10 b. All false reports are excluded, i.e., all zero-one 

variables in the second list whose subscripts have exactly 

one nonzero index. 

T 

c. All variables z 7 j for which (i. +1 , i. = (0, 0, 0) 

- L Jfc- L Jt+l'"- L Jk+N K 1 K N 

and Z{p) n T k (l k ) ={0) for p=k-K, k where K>0 is user 
15 specified. Thus the track T k {l k ) is terminated if it is not 

observed over K+J+l frames. 
Given the enumeration (4.15), one now fixes the assignments on 
the all index sets up to and including the (A:+l) index sets. 



20 T 11 ) 
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Thus one can now formulate the assignment problem for the 
next problem exactly as in (4.14) but with k replaced by k + 1. 
Thus, we do not repeat it here. 
5 The Initial Step. Here is one explanation for the initial 

step. First, assume that N=I+J. In this case, we start the 
track initiation with a solution of (3.5). Let 

{i 1 (i 2 ),i 2 (i 2 )...,i N (i 2 ),in 1 (i 2 ) Yi^o 

be an enumeration of the solution set of (3.5), i.e., those 
zero-one variables z. . . . . =1, including z 00 ... 0 =l 

1 ' 2 ' 2 ^ 2 ' " N+l ' 2' 

10 corresponding to I 2 =0, but excluding all those zero-one 
variables that are assigned to one and correspond to false 
reports (i.e., there is exactly one nonzero index in the 
subscript of zi^- ~i N+1 ) , all those zero-one variables that are 
assigned to one and correspond to tracks that initiate on 

15 frames higher than T+2 . Then we fix the data association 
decisions corresponding to the reports in our list of tracks 
prior to and including frame Jc+l = 1+2 . This defines the k for 
problem (4.4) and one can then continue the development by 
adding a frame to the window as in the general case. 

20 If J+J>N, then one possibility is to start the process 

with N+l frames, and assuming J^N, proceed as before replacing 
J by N-J for the moment, and continue to add frames without 

176 



NUMERI .01USC2 

lopping off the first frame in the window until reaches a 
window of length T+J+l. Then we proceed as in the previous 
paragraph . 

If I+J<N, then one can solve the track initiation problem 
5 (3.5), formulate the problem with the center of the window at 
Jc+l = N+l-J, enumerate the solutions as above, and lop off the 
first N-J-I frames. Then, we proceed just as in the case 
I+J=N. 

A primary objective in this work has been to demonstrate 

10 how multidimensional assignment problems arise in the tracking 
environment. The problem of track initiation and maintenance 
has been formulated within the framework of a moving window 
over the frames of data. The solution of these NP-hard, noisy, 
large scale, and sparse problems to the noise level in the 

15 problem is fundamental to superior track estimation and 
identification. Thus, one must utilize the special structure 
in the problems as well take advantage of special information 
that is available. Since these moving windows are overlapping, 
there are some algorithm efficiencies that can been identified 

20 and that take advantages of the overlap in the windows from one 
frame of reports to the next. Here is an example of the use of 
a primal solution of one problem to warm start the solution of 
the next problem in the sequence . 

Suppose we have solved problem (4.4) and have enumerated 

25 all those zero-one variables in the solution of (4.4) as in 
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(4.5). Add the zero index l k+1 =0 , so that the enumeration is 



L (5.1) 

, jt+i 

Vi =0 



With this enumeration one can define the cost by 



Q — C 

1 k+l^k+l+N ^Jfc^Jt+l' ' ^Jfc+l* ' '"' 1 N+l f 2 Jt+l+N 



(5.2) 



5 and the two dimensional assignment problem 

3> = Minimize 22 Y, c 7 ? z 7 ? =V 7 {z 2 ) 

Subject to: J. z 7 _i 7 = i r (5.3) 

■""k+l+N u 

Y z 2 ■ = 

1 =n "^Jt+l^Jt+l+W ^' ^Jt+l+N ^' ""'"jt+l+N' 
■"-k+l u 

2 



Let be an optimal or feasible solution to this two- 
dimensional assignment problem and define 
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(5.4) 



1 lf ^Jt+l'""' 1 k+N 




k+1' k+N 





and w, . =1 



jt+i 



or if d k+1 ,i k+1+N : 



0 otherwise. 



This need not satisfy the constraints in that there are 
usually many objects left unassigned. Thus, one can complete 
the assignment by using the zero-one variables in (4.4) with k 
5 replaced by k+1 with exactly one nonzero index corresponding to 
any unassigned object or data report. 

For the dual solutions, the multipliers arising from the 
solution of the two dimensional assignment problem (5.3) 



10 These are good initial values for use in a relaxation scheme 

[11,12] . Finally, note that one can also develop a warm start 

for problem (4.14) in a similar fashion. 

Multidimensional assignment problems govern the central 

problem of data association in multisensor and multitarget 
15 tracking, i.e., the problem of partitioning observations from 

multiple scans of the surveillance region and from either 



corresponding to the second variable, i.e 




179 



NUMERI .01 USC 2 

single or multiple sensors into tracks and false alarms. This 
fundamental problem can be stated as 



Minimize 



Z< • • • Z-r C U...i„ Z i,...t K 



i =0 1 " 1 

N u 



Subject to: 



2-/ • • ♦ 2^ 2~i ' • ' 2~i - • ' z i,...i =i 

for i p =l,...,M p and p=2,...,N-l, 
£ ••• S z i,-i„ = 1 ' V 1 '-'^' 



i 1 = 0 



z. . e {O, l} for all i, i 



(1.1) 



5 for i k =l,...,M k and k=2,...,N, 

zix-l^x e{0,l} for all i 1# ».,i N+1 
where c 0 ... 0 is arbitrarily defined to be zero and is included for 
notational convenience. One can modify this formulation to 
include multiassignments of one, some, or all of the actual 
10 reports. The zero index is used in representing missing data, 
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false alarms, initiating and terminating tracks. In these 
problems, we assume that all zero-one variables zi 2 -i w with 
precisely one nonzero index are free to be assigned and that 
the corresponding cost coefficients are well-defined. (This is 
5 a valid assumption in the tracking environment.) Although not 
required, these cost coefficients with exactly one nonzero 
index can be translated to zero by cost shifting without 
changing the optimal assignment. Finally, this formulation is 
of sufficient generality to include the symmetric problem and 

10 the asymmetric inequality problem. 

The data association problems in tracking that are 
formulated as (1.1) have several characteristics. They are 
normally sparse, the cost coefficients are noisy and the 
problem is NP-hard, but it must be solved in real-time. The 

15 only known methods for solving this NP-hard problem optimally 
are enumerative in nature, with branch -and -bound being the most 
efficient; however, such methods are much too slow for real- 
time applications. Thus one must resort to suboptimal 
approaches. Ideally, any such algorithm should solve the 

20 problem to within the noise level, assuming, of course, that 
one can measure this noise level in the physical problem and 
the mathematical method provides a way to decide if the 
criterion has been met. 

There are many algorithms that can be used to construct 

25 suboptimal solutions to NP-hard combinatorial optimization 
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problems. These include greedy (and its many variants), 
relaxation, simulated annealing, tabu search, genetic 
algorithms, and neural netork algorithms. For the three 
dimensional assignment problem, Pierskalla developed the tri- 
5 substitution method, which is a variant of the simplex method. 
Frieze and Yadegar introduced a method based on Lagrangian 
relaxation in which a feasible solution is recovered using 
information provided by the relaxed solution. Our choice of 
approaches is strongly influenced by the need to balance real- 

10 time performance and solution quality. Lagrangian relaxation 
based methods have been used successfully in prior tracking 
applications. An advantage of these methods is that they 
provide both an upper and lower bound on the optimal solution, 
which can then be used to measure the solution quality. These 

15 works extend the method of Frieze and Yadegar to the 
multidimensional case. 

Probabilistic Framework for Data Association. (ABP) 
The goal of this section is to explain the formulation of 
the data association problems that governs large classes of 

20 data association problems in centralized or hybrid centralized- 
sensor level multisensor/multitarget tracking. The 
presentation is brief; technical details are presented for both 
track initiation and maintenance in the work of Poore. The 
formulation is of sufficient generality to cover the MHT work 

25 of Reid, Blackman and Stein and modifications by Kurien to 
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include maneuvering targets. As suggested by Blackman, this 
formulation can also be modified to include target features 
(e.g., size and type) into the scoring function. The recent 
work of Poore and Drummond significantly extends the work of 
5 this section to new approaches for multiple sensor centralized 
tracking. Future work will involve extensions to track-to- 
track correlation. 

The data association problems for multisensor and 
multitarget tracking considered in this work are generally 
10 posed as that of maximizing the posterior probability of the 
surveillance region (given the data) according to 



where Z N represents N data sets, y is a partition of indices 
of the data (and thus induces a partition of the data) , r* is 

15 the finite collection of all such partitions, r is a discrete 
random element defined on r*, y° is a reference partition, and 
P(r=y|Z N ) is the posterior probability of a partition y being 
true given the data Z N . The term partition is defined below; 
however, this framework is currently sufficiently general to 

20 cover set packings and coverings. 




(2.1) 
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Consider N data sets Z (k) (k=I, ...,N) each of M k reports 




and let Z N denote the cumulative data set defined 



by 



Z(k) = 




iu=l 



and Z w = 



Z(l) ,...,Z(N) 



(2.2) 



respectively. In multisensor data fusion and multitarget 
5 tracking the data sets Z(k) may represent different classes of 
objects, and each data set can arise from different sensors. 
For track initiation the objects are measurements that must be 
partitioned into tracks and false alarms. In our formulation 
of track maintenance, which uses a moving window over time, one 

10 data set will be tracks and remaining data sets will be 
measurements which are assigned to existing tracks, as false 
measurements, or to initiating tracks. In sensor level 
tracking, the objects to be fused are tracks. In centralized 
fusion, the objects may all be measurements that represent 

15 targets or false reports, and the problem is to determine which 
measurements emanate from a common source . 



convenience in representing tracks, we add a zero index to each 



We specialize the problem to the case of set partitioning 



defined in the following way. 



First, 



for notational 
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of the index sets in (2.2), a dummy report z to each of the 
data sets Z(k) in (2.2), and define a "track of data" as 

1 N \ 

z . , z 1 where @@i k and z k can now assume the values of 0 and 

Zq, respectively. A partition of the data will refer to a 
5 collection of tracks of data wherein each report occurs exactly 
once in one of the tracks of data and such that all data is 
used up; the occurrence of dummy report is unrestricted. The 
dummy report z serves several purposes in the representation 
of missing data, false reports, initiating tracks, and 
10 terminating tracks. The reference partition is that in which 
all reports are declared to be false. 

Next under appropriate independence assumptions one can 

show 

Ptr-vlz") , = Ly= _ n (2 . 3) 
P(r = Y °|z N ) (i 1 -i M )eY ^ 

15 Li^.-iH is a likelihood ratio containing probabilities for 
detection, maneuvers, and termination as well as probability 
density functions for measurement errors, track initiation and 
termination. Then if ci 1 ...i N =-lnLi 1 ...i N , 
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(2.4) 



Using (2.3) and the zero-one variable zi^-i^l if (i x , 
Y and 0 otherwise, one can then write the problem (2.1) 
following N-dimensional assignment problem: 



..,i N ) e 
as the 



Minimi 



Subject to £ z =1 

E . . z. . =i 
y\ . z. . =i 



(i^l,...^) 

(i 2 -l,...,M 2 ) 



(2.5) 



(i p = l,...,M p andp = 2,...N-l) 
Eii...i z..... =1 (i N = l,...,M N ) 

z L t 6(0,1} for all i 1/ -,i w . 



where c 0 ... 0 is arbitrarily defined to be zero. Here, each group 
of sums in the constraints represents the fact that each-non- 
dummy report occurs exactly once in a "track of data." One can 
modify this formulation to include multiassignments of one, 
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some, or all the actual reports. The assignment problem (2.5) 
is changed accordingly. For example, if z k is to be assigned 
no more than, exactly, or no less than n k times, then the w = l" 
in the constraint (2.5) is changed to <, =, >n k , " respectively. 

5 Modifications for group tracking and multiresolution features 
of the surveillance region will be addressed in future work. 
In making these changes, one must pay careful attention to the 
independence assumptions that need not be valid in many 
applications . 

10 Expressions for the likelihood ratios Lii,...,^, can be 

found in the work of Poore. These expressions include the 
developments of Reid, Kurien, and Stein and Blackman. What's 
more, they are easily modified to include target features and 
to account for different sensor types. In track initiation, 

15 the N data sets all represent reports from N sensors, possibly 
all the same. For track maintenance, we use a sliding window 
of N data sets and one data set containing established tracks. 
The formulation is the same as above except that the dimension 
of the assignment problem is now N+l. 

2 0 Overview of the Lagrangian Relaxation Algorithm. (ABP) ] 

Having discussed the N-dimensional assignment problem [(1.1)] 
(3.1) , we now turn to a description of the Lagrangian 
relaxation algorithm. The algorithm will proceed iteratively 
for a loop k~\ f . . ., N-2 . At the completion, there remains 

25 one two-dimensional assignment problem that provides the last 
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step which yields an optimal (sometimes) or near-optimal 
solution to the original N-dimensional assignment problem. In 
step k of this loop (summarized in Section [41 IV. 3 ) one starts 
with the following (Itf-7c+ [1] 1) -dimensional assignment problem 
with one change in notation. [ ] If Jc=l, then the index notation 
1 2 and L x are to be replaced by i 2 and M lf respectively. 

... . . N-k+l W-Jt+l 

Minimize >, c 7 z, * 

1 *' 1 Jt+l — 

Subject To z, N ? +1 * =1, l t = l f . . . , L tf (4.7) 

EN-* + l =1 i =1 M 

1 k 1 k+2* * * 

for i = 1, . . . ,M p and p = Jr+2, . . . 

[(3.1)] 

To ensure that a feasible solution of Equation ( [3] 4 . [1] 7) 

always exists for a sparse problem, all variables with exactly 

one nonzero index (i.e., variables of the form z i~o.*o f° r 

. . . / L k and z 0 w I* + o„. 0 for i p =l, . . . , M. and p=7c+l [2] , . . . , N) are 
p p 

assumed free to be assigned with the corresponding cost 
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coefficients being well-defined. This assumption is valid in 
the tracking environment (A. B . Poore . Multidimensional 
assignments and multitarget tracking, partitioning data sets, 
I.J, cox, P. Hansen, and B. Julesz, editors, DIMACS Series in 
5 Discrete Mathematics and Theoretical Computer Science, American 
Mathematical Society, Providence, R.I.. 19:169-198, 1995; A.B. 
Poore. Multidimensional assignment formulation of data 
association problems arising from multitarget tracking and 
multisensor data fusion. Computational Optimization and 

10 Applications, 3:27-57, 1994) . 

[Subsection 1 Section IV. 3.1 presents some of the 
properties associated with the Lagrangian relaxation of 
([3]4. [1]7) based on relaxing the last {N-k) -sets of 
constraints to a two-dimensional one. [ Al Section IV. 3. 2 

15 describes a new approach to the problem of recovering a high 
quality feasible solution of the original (N-k+1) -dimensional 
problem given a feasible solution (optimal or suboptimal) of 
the relaxed two-dimensional problem is described hereinafter. 
A summary of the relaxation algorithm is given [hereinafter] in 

20 Section IV. 3. 3 , [following whichl and in Section IV .3. 4, we 
establish the maximization of the Lagrangian dual (an important 
aspect of the relaxation procedure) to be an unconstrained 
nonsmooth optimization problem and then present a method for 
computing the subgradients . 

25 
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IV. 3.1 The Lagrangian Relaxed Assignment Problem \ . 1 

The (N-k+1) -dimensional problem ([3]4.[1]7) has (N-k+1) 

sets of constraints. A (M p +1) -dimensional multiplier vector 

(i.e., [u p ]u? E [MMp+i]!-^) associated with the [p-th]p th 

constraint set will be denoted by 

u p ={Uq, uf, . . . / u^) T with Uq = 0 being fixed for each p = k + 

p 

2, ... ,.W and included for notational convenience only. Now, the 
(N-k+1) rdimensionall -dimensional assignment problem ( [3] 4. [1]7.) 
is relaxed to a two-dimensional assignment problem by 
incorporating {N-k-1) sets of constraints into the obj ective 
function via the Lagrangian. Although any (N-k-1) constraint 
sets can be relaxed, we choose the last (N-k-1) sets of 
constraints for convenience. The relaxed problem is 
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<Wi< u * +2 ' s Minimize 0 N _ k+1 ( z w -* +1 ; u k+2 f . . . u N ) - 

... . , N-Jt+1 w-Jt + 1 

Minimize J. c i i . i z i ± »± 



p=*+2 i p =0 



W-Jt+1 



Minimize 



p=*+2 i p =0 



N-Jfc + 1 



- £ t < 

p=k-2 i p =0 p 



Subject To /J = If i t =lr • • 

1 k+l X k+2 *" 

E + 1 = I i =1 M 



(4.8: 



[(3.2)] 

One of the major stqps in the algorithm is the naximi Taticn of @ N _ k+1 (u* +2 , . . . , u N ) 
with respect to the multipliers {u k+2 , . . . , u N ) . It turns out 
that ® N - k+1 is a concave, continuous, and piecewise affine 
5 function of the multipliers ( u k+2 , . . . , u N ) , so that the 
maximization of & N _ k+1 is a problem of nonsmooth optimization. 
Since many of these algorithms [ [???]] require a function value 
and a subgradient of & N _ k+1 at any required multiplier value 
(u* +2 , . .., u N ) , we address this problem in the next subsection. 
10 We note, however, that there are other ways to maximize & N _ k&1 
and the next subsection just addresses one such method. 



191 



NUMERI . 0 1 USC2 



IV. 3. 2 Properties of the Lagrancrian Relaxed Assignment 
Problem \ . 1 

For a function evaluation of & N - k+1 * we show that an 
5 optimal (or suboptimal) solution of this relaxed problem 
([3]4. [2] 8) can be constructed from that of a two-dimensional 
assignment problem. Then, the nonsmooth characteristics of 
®N-k+i are addressed, followed by a method for computing the 
function value and a subgradient. 



10 



Evaluation of <ft 



N-k+l 



Define for each {l k ,i k+1 ) an index 



<J* + 2/ »■/ Jn) =(jjt + 2(-Zjfc/ijt + i) /-/Jw(Ijt/ijt + i) )_and a new cost function 



■ L k ± k+l 



by_ 



argmin^ 



i + £ «f U = 0,1,...,M 



p=Jt+2 p 

and p= Jr+2, . . . ,N 



N-k+l 



p=k*2 



-00" 



Minimum < 



°' C 00i i + S "if' 



(4.9) 



[(3.3)] 
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Given an index pair {l k ,i k+1 ) , (j k+2 f - • - / j N ) need not be unique, 
resulting in the potential generation of several subgradients 
( [3] 4 . 1 [1] 7) . Then, 

£ ( U k+2 u N ) = 

Minimize <P N ^ k ^ 1 (z 2 ;u k+2 f . . . , u N ) 

Y, c?i *l l ~ Zl £ U F (4.10) 
Subject To 22 z ij^ =1 ' 1 k =1 ' • • - L k' 



As an aside, two observations are in order. The first is that 
5 the search procedure needed for the computation of the relaxed 
cost coefficients in ([3]4.[3]9) is the most computationally 
intensive part of the entire relaxation algorithm. The second 
is that a feasible solution z N ~ k+1 of a sparse problem 
([3]4. [1]7) yields a feasible solution z 2 of ([3] 4. [4] 10) via 
10 the construction 



± k 1 k+l 



1, if zf~**\ .. ti = 1 for some ( i , »• , i ) , 



0, otherwise. 

thus, there are generally solutions other than the one nonzero 
index solution. 

15 The following Theorem [3] 4.1 gives a method for evaluating 

^w-jt+i anc * states that an optimal solution of ( [3] 4. [2] 8) can be 
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computed from that of ( [3] 4 . [4] 10) . Furthermore, if the 
solution of either of these two problems is e -optimal, then so 
is the other. The converse is contained in Theorem [3] 4. 2. 

Theorem [3] 4.1. Let [c^Jwi be a feasible solution to the 
5 two^ dimensional assignment problem ( [3] 4. [4] 10_) and define [of' 



^w^r^Vw if(i A+2 , . . -(j k [(3.5)] 

and (1 } (4.11) 

and (1 } 

. rr N-k + 1 , V 

if C 00 W i N + L 



10 



N-Jfc+1 _ - 
W 00i -i ~ L 



15 



w-Jt+i A . ^ w-Jt+i t v 

"ooi^-i, = 0 2. (4.11) 

Then [ of' k+1 1 w"'** 1 is a feasible solution of the Lagrangian 

relaxed problem and 

20 If, in addition, [o 2 ] is optimal for the two [ ] ^dimensional 
problem, then [gf* k+1 7 w 1 *^* 1 is an optimal solution of the relaxed 

problem and <P N _ k+1 ( [u k+2 ] t££, [i^Jif) = <p [N _ 

lzk+1 ( [u k+2 ]i££, . . ., [V s ] if) . [ 
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Proof. Let rco 1 *' 1 " 1 ] w N ' k+1 and [<o 2 ]ud be as in the hypotheses, and 

let <p„. k+1 ( [co N - k+1 ] n£±i/ [if^Jv^l, [if] if.) 

and_ 

(z 2 ; i^^ltl, . . . , f j£) 

5 denote the objective function values of ([3]4.[2]8) and 
( [3] 4 . [4] 10) , respectively. Direct verification shows that 
rco N - k+1 1 w"-** 1 satisfies the constraints in ( [3]4. [2] 8) and 

<p N . k+1 ( [co N " k+1 ] [iS^lutl; [if ]ifl= <p [N . 

,^ k+1 ( [m 2 ] wi / f^Ju^,..., [if]if±._ 

10 For the remainder of the proof, assume that [co 2 ]wf. is optimal 
for ( [3] 4. [4] 10) . Let ry N ' k+1 l x w ' JctI satisfy the constraints in 
([3]4. [2] 8) and define 

*i 2 i =Ei i ^ W ? +1 i -i for (1^,1^)* (0,0), 
15 x 0 2 0 = l if" (l k ,i k+1 ) = (0,0) and 



(4.12) 



*' 2 w p 



x 0 2 0 = 0 if (1^1^)^(0,0) and 



20 



** 2 " p=Jc+2 p 



_Note that [x 2 ]2d satisfies the constraints in 



( [3] 4 . [4] 10) and 
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<P N . k+1 ( [x N " k+1 ] 2d±i; l^liJtl, • - - , 1^1 if.) > Sis- 

lzk+1 ( [X 2 ] 2d; [xf" 2 ]l£l [if)> N . k+1 (J; u* +2 xf) = <P„- k+1 (of- 

*_<V* +1 (w 2 : u*+ 2 u") 

5 = d> rIril (^- k+1 :u k+2 if) 

for any feasible solution r- x N " k * 1 1 x"'^ 1 of the constraints 
( [3] 4. [2] 8) . This implies rco N ' k+:L l w"~ k+1 is an optimal solution of 
( [3]4. [2] 8) . Next, 

<P N . k .! ( [u k+2 ] j£f, . . . , [if] if) = S N . k+1 

10 ( [u k+2 ] u^, . . . , [if]if)_ 

follows immediately from 

tfw ( [c>> N - k+1 ] is^±i; i^ 2 ]vtl, ., if] u?) = <p n . 

k+1 ( [O) 2 ]^/ [v** 3 ]^ , [f]f)_ 

for an optimal solution [co 2 ]^ of ( [3] 4. [4] .10) . [ 
15 ]_With the exception of one equality being converted to an 

inequality, the following theorem is a converse of Theorem 
[314.1. 

Theorem [3] 4.2. Let fof'^ 1 7 if'"* 1 be a feasible solution to 
problem ([3] 4. [2] 8) and define [cj 2 ]*? by 
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(4.13) 



"11 = E "?T\ -i for (l k ,i k .,)* (0,0) , 

1 k 1 k+l . . ■ L k :L k+l 1 k+2 x tf * * + l 

2 Jt+2 " * 2 M 

^ 0 2 0 = 1 if (l Jk ,i Jk+1 ) = (0,0) and 

c oo£-i H + £ u £ *° for some ^w™, i ») ' 

p=Jt+2 p 

^ 0 2 0 = 0 if (1^,1^) = (0,0) and 

Ci. + £ < >0 for aIi • 

* 2 " p=*+2 p 



10 Then [afjv? a feasible solution of the problem ([3] 4. [4] 10) and 
<P N . k+1 ( [co N " k+1 ] w^±f; fy* +2 Jj£!, . . . , fi/M) > 4- 

[co 2 ]ivf / [iS* 2 ]^, fi^Jif) ._ 

Jf, in addition, /" o?- k+1 ]v£±i is optimal for ([3] 4. [2] 8) , then 
[co 2 ]^ is an optimal solution of ( [3] 4. [4] 10) , _ 
15 0 w _ k+1 ( [w N - k+1 ] u£±i; /V +2 Ji£f, . . . , /"i/M) = 4. 

* + i ( [« 2 ] sd; [v** 2 ] and_ 

<P N . k+1 ([u k+2 ]i££ ,[iflif) = <p N . 



( [u k+2 ] u^f, . . . , fiA7u?) . /" 



Jc+2 
] 

20 Proof: Let rco N " k+1 1 w"'** 1 and [co 2 ].^ be as in _the_ hypotheses, 

and let <p„_ M ( [co N " k+1 ] v£±i; [if+ 2 l\£±, . . . , [i^Jif) and <p N . 

k+i ( [w 2 ]isd; fu** 2 ! u k+2 , . . . , [ifjif.) denote the objective function 
values of ([3]4.[2]8) and ( [3] 4 . [4] 10) , respectively. Direct 
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verification shows that [co 2 1 w 2 satisfies the constraints in 
( [3] 4. [4] 10) and 

<A™ ( [co N - k+1 ] vf^i; f^ +2 J i££, . . . , [if] if) > 0 S _ 

k +A [co 2 ]^; [if^Ju^ , [if] if) ._ 

5 For the remainder of the proof, assume that fo) N ' k+1 1 vf- k+1 is 
optimal for ([3]4.[2]8) and construct [co 2 ] wf as above. Using 
Theorem [3] 4.1, construct w N ' k+1 [N " k+1] by replacing rco N ' k+1 l vf~ k+1 in 
([3] 4. [5] 11) with w N ~ k+1 tN " k+13 . Then, from that theorem 
( " N - k+1 [N ~ k+1] ; [v" +2 ] ±i, if) -[] = [] <P*- k+ i ( [co 2 ] jsd/ [ 

10 if +2 ] j£f , []..., [if) < <P N . k+1 (of~ k+1 ; if* 2 , ...,if). ]xfl_ 

^. k+1 {^- k+1 ;u k+2 f . ..,u»). 

Optimality of [co N_k+1 1 vf~ k+1 then implies the last inequality is 
in fact an equality. This proves the theorem. 

The Nonsmooth Optimization Problem. Next, we address the 

15 nonsmooth properties of the function @ N ^ k+1 { l"u k+2 1 u k+2 , . . . , [if] if) 
as explained in the following theorem. 

[Theorem 3.4 ] Theorem 4,3, (G.L. Nemhauser and L.A. 
Wolsey. Integer and Combinatorial Optimization, Section II. 3.6. 
Wiley, New York, NY, 1988) . Let <P N . k+1 be as defined in 

20 ( [3] 4. [2] 7) , let [N N _ k+1 ] V N _ k+1 (z N ' k+1 ) be the [object] objective 
function value of the (N-k+1) -dimensional assignment problem in 
equation ( [3] 4. [1] 7) corresponding to a feasible solution z N ' k+1 
of the constraints in ([3]4.[1]Z), and let 
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z N - k+1 be an optimal solution of ([3]4.[1]7). Then, <P N _ 
kil ( [u** 2 ] ^* 2 , . . . , [jf]u N ) is piecewise affine, concave and 
continuous in { \ u k+2 1 u k + 2 , . . . , [xf*]\F} and 

5 0 (u k+2 u N )<V (z N ~ k+1 )<V lz N ~ k+1 ) ^ 4 ' 

Most of the algorithms for nonsmooth optimization are based on 
generalized gradients, called subgradients , given by the 
following definition for a concave function. 

10 

Definition [3]4.[5]1. At u = ( u* +2 , . . - , u N ) the set d& N _ k+1 (u) is 
called a suhdifferential of & N „ k+1 and is defined for the 
concave function & N _ k+1 (u) by 

15 d<P N _ k+1 (u)={geR M ^ +1 x ... x IR^ 1 1 $ N _ k+1 (w) - <£ w _* +1 (u) * (4 ' 15) 
g T (w-u) for all we 1 w * +2+1 x...xR m » +1 }, 

inhere 

g^ = fc^=u;f = 0 are all permanently fixed. (Recall that these were 
used for notational conveience only.) A vector ge d& N _ k+1 is 
20 called a subgradient. 

There is a large literature on such problems, e.g. , (J. - 
B. Hiriart-Urruty and C. Lemarechal . Concex Analysis and 
Minimization Algorithms I& II. Springer- Verlag, Berlin, 1993; 
25 K.C. Kiwiel. Methods of descent for nondif f erentiable 
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optimization. In Lecture Notes in Mathematics 1133. A. Dold and 
B. Eckmann, eds. Springer-Verlag, Berlin, 1985; C. Lemarechal 
and R. Mifflin. Nonsmooth Optimization. Pergamon Press, Oxford, 
UK, 1978; B.T. Polyak. Subgradient method: 
5 A survey of Soviet research. In C. Lemarechal and R. Mifflin, 
editors, Nonsmooth Optimization, pages 5-29, N.Y. , 1978. 
Pergamon Press. ; H.Schramm and J. Zowe. A version of the bundle 
idea for minimizing a nonsmooth function: Conceptual idea, 
convergence analysis, numerical results. SI AM Journal on 

10 Optimization, 2, No. 1 : 121-152 , February, 1992; N.Z. Shor. 
Minimization Methods for Non-Pi fferentiable Functions. 
Springer- Verlag, Nem York, 1985; P.Wolfe. A method of conjugate 
subgradients for minimizing nondif f erentiable functions. 
Mathematical Programming Study. 3:145-173, 1975), and we have 

15 tried a variety of methods including subgradient methods (N.Z. 
Shor. Minimization Methods for Non-Pi fferentiable Functions. 
Springer-Verlag, New York, 1985) and bundle methods ( J. -B. 
Hiriart-Urruty and C. Lemarechal. Convex Analysis and 
Minimization Algorithms I& II. Springer-Verlag, Berlin, 1993; 

20 K.C. Kiwiel. Methods of descent for nondif f erentiable 
optimization. In Lecture Notes in Mathematics 1133. A. Dold and 
B . Eckmann , eds. Springer-Verlag, Berlin, 1985) . Of these, we 
have determined that for a fixed number of nonsmooth iterations 
(e.g., twenty), the bundle trust method of Schramm and Zowe 

25 (H.Schramm and J. Zowe. A version of the bundle idea for 
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minimizing a nonsmooth function: Conceptual idea, convergence 
analysis, numerical results. SI AM Journal on Optimization, 2. 
No. 1:121-152, February, 1992) provides excellent quality 
solutions with the fewest number of function and subgradient 
5 evaluations, and is therefore our currently recommended method. 
An Algorithm for Computing <2> w _ Jt+1 and a Subgradient 
Most current software for maximizing the concave function 
^w-jt+i re< ? u i res the value of the function and a subgradient at 
a point (u k * 2 , . . . , u N ) . Based on the previous two sections, 
10 this can be summarized as follows. 

1. Starting with problem ([3]4.[1]7), form the relaxed 
problem ( [3] 4 . [2] 8) ; 

2. To solve ([3]4. [2] 8) optimally, defined the two^ 
dimensional assignment problem ([3]4. [4] 10) via the 

15 transformation ( [3] 4. [3] 9) ; 

3. Solve the two-dimensional problem ([3] 4. [4] 10) 
optimally; 

4. Reconstruct the optimal solution, say ft N ~ k+1 , of 
([3] 4. [2] 8) via equation ([3] 4. [5] 9) as in Theorem [3] 4.1; 

20 5 

V Compute 
1 k 1 k+\'" 1 N the 



♦ e. £»d.. e •ir.i.-i 



25 



p=*+2 i p =0 p 



function 



201 



[(3.10)] 
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6 . A subgradient is given by 



dul 



for i p = 1 / • • • / M p and p = k+2, . 



(4,17) 



[ (3.11)] 

Several remarks are in order. First, the optimal solution of 
the minimization problem ( [3] 4. [2] 8.) is required before one can 
remove the minimum sign, replace z N ~ k+1 by the minimizer i^ w " jc+1 

10 and differentiate with respect to to obtain a subgradient 
as in ( [3] 4 . 1 [1] 7) . If ti**~ k+1 is an approximate solution of 
( [3] 4. [2] 8) , then the subgradient and function values are only 
approximate with [accruacy] accuracy depending on that of ^~ k+1 . 
Although one can evaluate the sums in ([3]4.1[0]6) and 

15 ([3]4.1[1]7) in a \ straightforward! straight forward manner, 
another method is based on the following observation. Given 
the multiplier vector {u k+2 , . . . , u N ) , let 

{UjtUjt + i> ' ) } Ljc+1 be an enumeration of indices {l k ,i k+1 | of 

w 2 (or the first 2 indices of w N ' k+1 constructed in equation 
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[3 ] J_4_. [ 5 ] 9 ) ) such that wf , , )± (1 v=l for 

and (l k (l kil ) . = (0.0) for l All =0 regardless of whether 

Wq 0 =1 or not. (The latter can only improve the recovered 
5 feasible solution.) The evaluation of the bracketed quantity 

in ([3]4.1[1]7) for a specific i p * 1 and rp=k+2l p = k+2 N 

is one minus the number of times the index value i p appears in 
the ( fp-k+1 7 p-k+1 ) th position of the ( fN-k+ll N-k+1 ) (tuple) -tuple 
in the list {l k d k+1 ), i w (i w )J W M..,j| w • 

10 Finally, we have presented a method for computing one 

subgradient . If ij**'** 1 is the unique optimal solution of 
( [3] 4. [2] 8.) , then ® N . k+1 (u) is dif f erentiable at u, g is just the 

gradient at u, and the subdif f erential d® N - k+1 ( u ) = { 7) is a 

singleton. If, corresponding to u,_^~ k+1 is not unique, then 

15 there are finitely many such solutions, say 
(1) , . . . , (m) } with respective subgradient s 

{g(D , . . . , g(m) } , [ ] the subdif f erential d<P{u) is the convex hull 
of {g(l) , . . . ,g(m) } (J.L. Goffin. On convergence rates of 
subgradient optimization methods, mathematical programming. 

20 Mathematical Programming, 13:329-347, 1977) . These nonunique 
solutions arise in two ways. First, if the optimal solution of 
the two[ ] ^dimensional assignment problem is not unique, then 
one can genera [l]te all optimal solutions of the two [ ]^ 
dimensional assignment problem ([3] 4. [4] 10) . Corresponding to 

25 each of these solutions and to each index pair (l k , i k+1 ) in each 
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solution, compute the indices ( j k+2 * • • • 9 j N ) i n ( [3] 4 . [3]9) . If 
these j's are not unique, then we can enumerate all the 
possible optimal solutions tf*'** 1 of ([3]4. [2]8). Given these 
solutions, one can then compute the corresponding subgradients 
5 from ( [3] 4.1 [1] 7) . [ 

]_In most nonsmooth optimization algorithms, one uses an 
[epsilon-subdif f erential] e-subdif f erential . 

Definition [3]4.[x]2. At u = (u k+2 , . . . , if) the set d e $ N _ k+1 {u) is 
10 called an [epsilon subdiff erential] e-subdif f erential of ^ N _ k+1 
and is defined for the concave function * w _ Jt+1 (u) f (u) ] by 

3 e <V* + i ( «) = i g e M M ^ +1 x ... x r* +1 1 0 N _ k+i ( w) - 0 N _ k+1 ( U ) * 
g T (w-u) + 6 for all w el w *- 2+1 x . . . * , 

where g^=^=u^ = 0 are all permanently fixed. (Recall that 
these were used for notational convenicence only.) A vector 
g e dgd^.^j (u) is called [ a subaradientl an e-subcrradient . 

15 

[HOW DO YOU COMPUTE g 6 d<J> N _ k+1 (u) ] 

IV. 3. 3 The Recovery Procedure. 
The next objective is to explain a recovery procedure, i.e., 
20 given a feasible (optimal or suboptimal) solution [co 2 ]wf of 
(C3]4. [4] 10) (or [co N - k+1 ] ss^LUll of ([3]4.[2]8) constructed via 
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Theorem [3]4.1), generate a feasible solution z^' k+1 of equation 
( [3] 4. [1] 7) which is close to [to 2 ] wf in a sense to be specified. 
We first assume that no variables in ( [3] 4. [1] 7) are 
pre [s] assigned to zero; this assumption will removed shortly. 
5 The difficulty with the solution [co N " k+1 ] in Theorem [3] 4.1 

is that it need not satisfy the last ( fN-k-1 W-k-1 ) sets of 
constraints in ([3]4. [1] 7) . (Note, however, that if [co 2 ]^ is 
an optimal solution ([3]4. [4] 10.) and rco N " k+1 1 as constructed 

in Theorem [3] 4.1, satisfies the relaxed constraints, then [to N ~ 

10 k+1 1 w"-^ 1 is optimal for ( [3] 4 . [1] 7) . ) The recovery 

fproceudrel procedure described here is designed to preserve the 
0-1 character of the solution [co 2 ]wf of ([3]4. [4] 10) as far as 
possible: If [] [o>] i [ =1 and ] =1 and l^*[0or]0 or 1^*0, 
the corresponding feasible solution z* _jc+1 of ([3]4.[1]7) is 

15 construed so that z 1 w "/ +1 „ i =l_for some (i k+2 , . . . , i N ) . By this 

k k+1 N 

reasoning, variables of the form z^"^..^ can be assigned a 
value of one in the recovery problem only if [o>jWq 0 = 1. 
However, variables z^"^..^ will be treated differently in the 
recovery procedure in that they can be assigned either zero or 
20 one independent of the value [of] [co^o- This increases the 

feasible set of the recovery problem, leading to a 

potentially better solution. 

[ ]Let { d k (l k+1 ) , i k+1 d k+l ) Y k+1 _be_ an enumeration _of 
indices {l k/ i k+1 } of [co 2 ]wf (or the first 2 indices of [co N ~ 
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k+1 1 u^ -jc+1 constructed in equation ([5)) such that for and for 
l k+1 =0] 4.11) ) such that 

w i a ) i a ) =1 for 

regardless of whether [co^^l or not. To define the ±N-k)_z 
dimensional assignment problem that [resores] restores 
feasibility, let 

r w_1 = c N for k=l 

N-* _ N-k+1 fA p v 

c l i i ~ c l (1 )i (1 )i i (4.10? 

= c w • .... 

for k*2 and 1 = 0, . . . , L , 



where 



10 for m = 2, . . .Jc+1 and where [o]^l denotes function composition. 
For example^ 1 23 = 1 2 {1 3 ) and 1 234 [m } [l 2 o 1 3 (1 4 )]= 1 2 [ (]£l 3 
1 2 ±1 3 11 4 ) ) . 

Then the {N-k) -dimensional assignment problem that 
restores feasibility is [ (3.14)] 
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kjt> • - N-k N-k 

Minimize ), ^ 7 z 7 

, . Jt+l Jr+2 * * ' N - L lc+l- L Jt+2" " • X N 
J, * + l 1 *+2' ' " 2 N 

Subject To £ z w M ...i„ =1 ' 1 m = 1 '""V 

^Jk+2 ' • ' 2 W 

EN"* =1 7 =1 M 



1 k 1 k*l ■ • • i p-l i p+l • • • 1 N 



N-k 



= 1 



for i=l, . . .;M and p-k+3, . . . , N-l, 



L 2 Jt+2* ' • 2 W-1 



(4.20) 



Let Y be an optimal or feasible solution to this [N - kl (N-k) - 
dimensional assignment problem. The recovered feasible 
solution 2^ is defined by [ (3.14) ] 



z ± i i i = ^ 

dL l" L 2 x 3 X W 



1, if i^i^lz k+1 ) , i 2 = i ? (-i2... Jk+1 )f 

i Jfl = i Jk*l^ I Jk + l> and Y l i i =1 ' 

0, otherwise. 



(4.21) 



Said in a different way, the recovered [reasible] feasible 
5 solution [z N ]z^ corresponding to the multiplier set 
{ ru k+2 l u k * 2 . is then defined by 
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1, if Y x . .... =1, 
Of otherwise, 



where l m .. k+1 is defined in ([3]4.1[3]9) and [o]«o^_ denotes 
function composition. 

The next objective is to remove the assumption that all 
5 cost coefficients are defined and all zero-one variables are 
free to be assigned. We first note that the above construction 
of a reduced dimension assignment problem ( [3] 4. [14] 20) is 
valid as long as all cost coefficients c^"* are defined and all 
zero-one variables in are free to be assigned. 

10 Modifications are necessary for sparse problems. If the cost 



coefficient c 7 ", / w w 7 - is well-defined and the zero-one 



variable z, w A +1 w , 7 w is free to be assigned to zero or 

one, then define cf~ k ± ± = ti w ...i as i n ([3] 4.1 [2] 8) 

with zT k ; i being free to be assigned. Otherwise, zT k i , 

15 is preassigned to zero or the corresponding arc is not allowed 
in the feasible set of arcs . To ensure that a feasible 
solution exists, we now only need ensure that the variables 

z i~*o-.q are free to be assigned for [l hll =0l l kll =0 ,1 L kil with 

the corresponding cost coefficients being well-defined. 

20 (Recall that all variables of the form zf^o and z^-o are 
free (to be assigned) and all coefficients of the form c^,*^ 1 
and 
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c o~oi)>-.o are well-defined for [l k =0] l k =0_, L k and [i p =0] ± P =Q , M p 
for p=k+l, ...,N. ) This is accomplished as follows: If the cost 
coefficient c lfc(WWW(M) is well-defined and Zj 4 ,i w ,i w ,j w ,w 
is free, then define c£"* 0 „ 0 = c^" { ^ l)ijt+l ( ijt+l )o.--o with z w>-o being 
free. Otherwise, since all variables of the form z^qSq and 
z oi* + J>-o are f r ^e to be assigned with the corresponding costs 
being well-defined, set c Ijk+i0 .„ 0 = c Ijt(ljM)00 _ 0 + c oijt+i(ijt+i)0 .„ 0 ^ where 
the first term is omitted if l k (l k+1 )=0 and the second, if 
ijt+i^+i^ 0 - For ^^Jt+i^ 0 and ijt+i^Jt+i^ 0 ' define c 0 w ;<f = Cq^' 1 . 



IV. 3. 4. The last step k = N-2. 

The description of the algorithm ends with rk=N+2l k = N+2 . 
The resulting recovery problem defined in ( [3] 4 . 1 [4] 0.) with 
[k=N-2 I Jc = N-2 is 



0 o = Minimize V] JZ cf L z\ ± =V„(z 2 



Subject To £ z i„,i„ =1 ' Vi = 1 ""'Vi' (4.22) 

i M =o "- 1 



15 Let Y be an optimal or feasible solution to this [2- 
dimensionall two-dimensional assignment problem. The recovered 
feasible solution zF is defined by 
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• L l- 1 2 2 3"'- 1 N 



1, if i x = i x {l 2 -vh^ ' L 2 = i 2 ^z-vh-d i i 3 = i 3 ( ly-v-lHUMERl .01USC2 
V2 = ^-2 < 2 W -2N-1> ' Vl = Vl ( Vl> Snd ^i, = 1 

or if- (Vi'V = <°'°) ' 
0, otherwise. 



( f3. 181 4. 23 ) 

Said in a different way, the recovered feasible solution is 
5 then defined by 



15 



20 



1, if y 2 . =1, (4.24) 

0, otherwise. 



[ (3.19) ] 

where l m .. N -! is defined in ([3] 4.1 [3] 9) and M "°" denotes 
function composition. To complete the description, let 
l(l N _ 1 (l N ),i N (l N )\ LN be an enumeration of indices {l^, i N \ of Y 
10 such that r Vi(Vf . M(V =l for (l^ (1 N ) , i w (2 N )) # (0, 0) and 
(V^V ,i M U N )) = (0,0) for 1 N =0 regardless of whether 3T 00 =1 or 
not. Then the recovered feasible solution can be written as 



1 f if i 1 = i 1 ( 1 2 J ,i 2 = i 2 ( l 2 „. w ) , i 3 = i 3 ( l 3 .J , ... , 
0, otherwise. 



(4.25) 



IV. 3 . 5 . The Upper and Lower [Bound. 1 Bounds 

The upper bound on the feasible solution is given by 



( [3]4.2[1]6) 



and the lower by [<I> W ] $ N (u 3 , u N ) for any multiplier value 
(u 3 ,..., fir^u^) . In particular, we have 
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0 N (u\...,u N ) I 0 N (u\...,u N ) Z V N (z N )z V N (2 N ), ([3]4.2[2]7) 

where (u 3 ,...,^) is any multiplier value, (u 3 ,...,^) is a 
5 maximizer of [® N ]0 N {u 3 ,...,u N ), is an optimal solution of the 
multidimensional assignment problem ([3]4. [1]7) and is the 
recovered feasible solution defined by ([3]4.2[0]5) . Finally, 
we conclude with the observation that V N {z) =V 2 (Y) where Y is the 
optimal solution of ( \3 . 171 4 . 23 ) used in the construction of z 
10 inequations ( [3] 4 . [18] 24) - ( [3] 4 . 2 [0] 6) . 

M IV.3.6. Reuse of Multipliers f . 1 

Since the most computationally expensive part of the 
algorithm occurs in the maximization of ^ N .^ +1 , the development 

15 of algorithms that can make use of hot starts or multipliers 
close to the optimal are fundamentally important for real-time 
speed. The purpose of this section is to demonstrate that the 
multiplier set obtained at stage kzl provide good starting 
values for those obtained at step Jc+1. 

20 Theorem [3]4.[xx]4. Let (u 3 ,...,u N ) denote a multiplier set 
obtained in the f maximum 1 maximi z at ion of & N [ (] (u 3 , u^) . Then 
this multiplier set satisfies the string of inequalities 
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& N (u 3 , ...,u N ) ^0 N . 2 (u 4 , . . .,u w ) s . . . z0 3 (u N ) Z0 2 ZV N (¥) , (4.28) 



[ . J 

where after the first step in the maximization of [<I> N ]@ N the 
multipliers are not changed in the remaining steps. 
Furthermore, to improve this inequiality, let 
(u 2+k ' N - k+1 ,...,u N ' N - k+1 ) denote a maximizer of & N _ k+1 ( u k+2 , u N ) . 
Then, we have 

^^(u 3 **'^** 1 , ...,u w ' w -* +1 ) 
z$ N _ k (u 3+k > N - k , ...,u N ' w "*) 

<L . . . *<P 3 (U N ' 3 ) Z$ 2 <L V N (f) . 



[ ( 3 . xxx) ] 

[ Inequalities depend on solution of 2d assignment problem.] 
Proof: Suppose we have a value of the multipliers (u k+2 , u N ) 
obtained in the maximization of 



<V* + i < u * +2 ' . • • , u») = Minimize <p N _ k+1 (z N ~ k+1 ; u k+2 , ...,u N ) (4.30) 
Subject To £ , Cw-i/ 1 ' 

E w-*+l =1 7 =1 M 



where [(3.xx) ] 
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^ t - N-k+1 . ,,/c+2 „ W\ 

W-Jt+l z N-k+l 



♦ £fc 

p=k+2 i p =0 



^k x k+l ' " * ^-l^p+l ' * • 2 N 



Z I i i " 1 



U i/ u i„ ^^^...i^ u i 



p=k + 2 



p=k+2 i p =0 



(4.31) 



These need not be the maximizer; however, we do assume that we 
have solved the above minimization problem optimally to 
evaluate @ N _ k+1 (u k+2 , u N ) . Just as in the definition of the 

[earlier ] recovery proble m discussed earlier , let 
5 {(l k (l k+1 ) ,i k+1 U k+1 ) )} L £\=o be an enumeration of indices 

{ l k , i k+1 [ ) ] J_ of the optimal solution w 2 of problem ([3]4. [4] 10) 

(or the first 2 indices of the solution fw N=k+1 1 v/*"** 1 of the 
relaxed problem ([3] 4. [2] 8) constructed in equation 

([3] 4_. [5] 11)) such that w iWi>*WW =1 for 

regardless of whether w^ 0 = l or not. 

Then, the final value of the objective function can be 
expressed as 
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Minimize <P N . k+1 {z N ~ k+1 ; u k+2 r u N ) (4.32) 
Subject to £ ZiX^^i^iwi, 5 ^ ^'-'Vi' 



where 

0 n ^ +1 (z N -* +1 ;u* +2 ,...,u N ) = (4.33) f(3.xxx)1 



- £ £ < . 



p=k+2 i p =0 



If now another constrain set, say the (Jc+2) set, is added to 
5 this problem, one has 
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0 N=k+1 (u k + 2 ,...u N ) = Minimize 0^ (z N "* +1 ; u* +2 , u N ) = (4. 34) 

Minimize £ . | ^www-i.* £ < f<<Wi*.i<W-i» 

- £ t < 

p=k+2 i p =0 p 

Subject tO £ <u^ 1 )i ttl (^ 1 ) W .i N =lr Vi = 1 '"' £ M' 

1 k+2' 1 N 

Ez N ~* +1 =1 i =1 AT 



Since the constraint £ ^UhIIhIWV 1 for 

I Jt+i :L Jt+3*" i N 

1^2=1 / ■••/ M k+2 is now imposed, the feasible region is smaller, so 
that one has 
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s <V k+1 (u* +3 ,...,u N ), 



where the fact that @ N = k+1 does not depend on the multiplier set 
u k+2 due ttie ^^3^ constraint set. Also, the last equivalence 
follows from the fact that & N = k+1 (u** 3 , [uFjy^) is the 

5 relaxed problem [ (3.xx)] (with 7c replaced by Jc+l) for the 
recovery problem [ (3.xx)] in step k of the above algorithm. 
Continuing this way, one can compute {u 3 ,...,uP) at the first step 
(7c=l) , fix them thereafter, and perform no optimization at the 
subsequent steps to obtain 
10 tf w (u 3 ,...,u N ) £* w _ 1 {uV..,u w ) *.»s* 3 <u w ) <0 2 =V N (2), 

where & 2 is defined in ([3]4. [17] 22) . 

To explain how to improve this inequality, let 
(u 2+ *' w ~* +1 ,...,u M ' w "* +1 ) denote a maximizer of ^_ k+1 ( u k+2 , u w ) . 
Then by the same reasoning one has I" ( 3 . xx) 1 the result as stated 
15 in the theorem. [Q.E.D.] 



Summary of the Lagrangian Relaxation Problem. (APB) 
V. SUMMARY OF THE LAGRANGIAN RELAXATION PROBLEM 
Given the multidimensional assignment problem 
20 [ (4.1)] 
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Minimize V, c. . z. 
i x — j. n 

Subject To Y\ z. . =1 (i=l,...M.), (5.1) 
i 2 i 3 ...i. 

V 3 ...i N 

X) z. . =1 (i=l,...M and p = 2, . . . , W-l) , 

. z i 1 ...i= 1 ' U N =1,..-M N ), 

z. e {0,1} for all i., .. . i, 

where c*. 0 is arbitrarily defined to be zero and is included for 
notational convenience and where the superscript N on both c 
and z is not an exponent, but denotes the dimension of the 
subscripts and the assignment problem as stated in ([2]3.5). 
5 In relaxed and recovery problems c 0 w . 0 need not be zero! In this 
problem, we assume that all zero-one variables z i^i N with 
precisely one nonzero index are free to be assigned and that 
the corresponding cost coefficients are well-defined. (This is 
a valid assumption in the tracking environment (A.B. Poore. 

10 Multidimensional assignments and multitarget tracking, 
partitioning data sets, I.J. Cox, P. Hansen, and B. Julesz, 
editors, DIMACS Series in Discrete Mathematics and Theoretical 
Computer Science, American Mathematical Society, Providence, 
R.I.. 19:169-198, 1995; A.B. Poore. Multidimensional assignment 

15 formulation of data association problems arising from 
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multitarget tracking and multisensor data fusion. Computational 
Optimization and Applications, 3:27-57, 1994) .) Although not 
required, these cost coefficients with exactly one nonzero 
index can be translated to zero by cost shifting (A.B. Poore 
5 and N. Riiavec. A lagrangian relaxation algorithm for 
multidimensional assignment problems arising from multitarget 
tracking. SI AM Journal of Optimization. 3, No. 3:544-563, 1993) 
without changing the optimal assignment. 

Having explained many of the relaxation features, it is 
10 now appropriate to present the Lagrangian relaxation algorithm, 
which iteratively constructs a feasible solution to the N- 
dimensional assignment problem ([4] 5.1). 

Algorithm [4] 5.1 (Lagrangian Relaxation Algorithm) [.] To 
construct a high quality feasible solution, [eno ted] denoted by 
15 z®, of the assignment problem ([4] 5.1), proceed as follows: 

0. Initialize the multipliers (u k+2 ,.~, u N ) , e.g., 

(u* +2 ,...,u w ) = (0,...,0) . 

For Jc-1, ...,N-2, do 

1. For the Lagrangian relaxed problem {[3] 4. [2] 8) from the 
20 problem ( [3]4. [1] 7) by relaxing the last (N-k-1) sets of 

constraints . 

2. Use a nonsmooth optimization technique to solve 

[ (4.2)] 
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Maximize {0 N _ k+1 \u p eM Mp+1 for p = fc+2,...,N (5.2) 

with Uq-0 being fixed], 



where & N _ k+1 ( u k+2 , u N ) is defined by equation ([3]4. [2] 8) . 
The algorithm [following problem d in Section IV. 3. [9)] 2 
provides one method for computing a function value and a 
subgradient out of the subdifferential at ( u k+2 , u N ) , as 
5 required in most nonsmooth optimization techniques. 

3. Given an approximate or optimal maximizer of ( [4] 5.2) , say 
{Q k+2 , ■», Q N ) y let [w 2 1 ft 2 denote the optimal solution of the 
two-dimensional assignment problem ( [31 4 . [4] 10) 
corresponding to this maximizer of & N _ k+1 ( u k+2 f u w ) . 
10 4. Formulate the recovery (N-k) -dimensional problem 
( [3] 4. 1 [3] 9) , modified as discussed at the end of 
subsection [ (1 IV. 3.3 f) 1 for sparse problems. At this 
stage, z®^ as defined in ( [3] 4. [15] 21) ^_ contains the 
alignment of the indices {i lf i k+1 } . 
15 Formulate the final two-dimensional problem ( [3] 4. [17] 22) and 
complete the final recovered solution z N as in ( [3] 4 . 2 [0] 5) • 
To explain the lower and upper bounds, let & N be as 
defined in ([3]4.[2]8) with k=l, let V N (z^) be the objective 
function value of the ^-dimensional assignment problem in 
20 equation ([2] 3.-5) corresponding to a feasible solution z 3 * of 
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the constraints in ([2] 3. 5), and let ^' k ^ b e an optimal 
solution of ([2] 3. 5). Then, 

<P N (u 3 ,...,u N ) <> V N {z N ) z V N (2 N ) ([4]5.3) 

is the desired inequality. 
COMMENTS : 

1. Step 2 is the computational part of the algorithm. 
Evaluating d^.jt+i and computing a subgradient used in the 
search procedure requires 99% of the computing time in the 
algorithm. This part uses two [ ] ^dimensional assignment 
algorithms, a search over a large number of indices, and 
a nonsmooth optimization algorithm. It is the second part 
(the search) that consumes 99% of the computational time 
and this is almost entirely parallelizable . 

2. In track maintenance, the warm starts for the initial 
multipliers for the first step are available. For the 
relaxation procedure, initial multipliers are available in 
step 2 from the prior step in the algorithm. 

3. There are many variations on the above algorithm. For 
example, one can compute a good solution on the first pass 
(7c=l) and not perform the optimization in (2) thereafter. 
This yields a great solution. Thus one can continue the 
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optimization at the first pass, and immediately compute 
quality feasible solutions to the problem. 
[ Lagrangian Relaxation Algorithm for the 3 -Dimensional 
Algorithm. (ABP) ] VI. LAGRANGIAN RELAXATION ALGORITHM 

5 FOR THE THREE-DIMENSIONAL ALGORITHM 

_ In this section, we illustrate the relaxation algorithm for 
the three [ ] ^dimensional assignment problem. Having discussed 
the general relaxation scheme, 



Minimize £ £ 2 c m z i±i 
Vo vo ivo 

Subject To £ £ z... =1, ^ = 1,...,*^, 

i 2 = 0 i 3 =0 123 

2 1 = 0 2 3 =0 1 2 3 

i 1= o i 2 =o ili2is J J 



( [5]6.1) 



z V 2 i 3 e(0,l) for all 1,1223. 
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To ensure that a feasible solution of ([5] 6.1) always exists 
for a sparse problem, all variables with exactly one nonzero 
index (i.e. , variables of the form z ± 00 , z Qi 0 , and z 00i ^ for 
i p -l f ...,M p and p-1,2,3 are assumed free to be assigned with the 
5 corresponding cost coefficients being well-defined]_. This 
assumption is valid in the tracking environment [ [**refPoa, 
**refPob] . 

] (A. B . Poore. Multidimensional assignments and multitarget 
tracking, partitioning data sets, I.J. Cox, P. Hansen, and B. 

10 Julesz, editors, DIMACS Series in Discrete Mathematics and 
Theoretical Computer Science, American Mathematical Society, 
Providence, R.I., 19:169-198, 1995; A. B . Poore. 
Multidimensional assignment formulation of data association 
problems arising from multitarget tracking and multisensor data 

15 fusion. Computational Optimization and Applications , 3:27-57, 
1994) . 

VI. 1. The Lagrangian Relaxed Assignment Problem f . 1 
The three [ ] ^dimensional assignment problem ([5]j6.1) has 

2 0 three sets of constraints and one can describe the relaxation 
by relaxing any of the three sets of constraints, the 
description here is based on relaxing the last set of 
constraints. A (M 3 +l) -dimensional multiplier vector 

(i.e., 3 m>M,+i) associated with the 3 tth3 — constraint set will 

25 be denoted by u 3 = |u 0 3 , u 3 , u^j T with u 0 3 = 0 being fixed for 
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notational convenience. (The zero multiplier u 0 3 =0 is used for 



notational convenience.) Now, the three-dimensional assignment 
problem ([5] 6.1) is relaxed to a two-dimensional assignment 
problem by incorporating last set of constraints into the 
5 objective function via the Lagrangian. Although any constraint 
set can be relaxed, we choose the last set of constraints for 
convenience. The relaxed problem is 




( t5] 6.2) 



<P 3 {u 3 ) = Minimize <p 3 (z 3 ;u 3 ) 



= Minimize 



±±± k 



+ U 




i 1= 0 i 2 =0 i 3 -0 



Subject to 




One of the major steps in the algorithm is the 
10 maximization of 0 3 ([u 3 ]ui) with respect to the multiplier 
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vector u 3 . It turns out that <P 2 is a concave, continuous, and 
piecewise affine function of the multiplier vector u 3 , so that 
the maximization of <P 2 is a problem of nonsmooth optimization. 
Since many of these algorithms require a function value and a 
5 subgradient of ® 3 at any required multiplier value (u 3 ) , we 
address this problem in the next subsection. We note, however, 
that there are other ways to maximize ^ 3 and the next subsection 
addresses but one such method. 

10 VI. 2. Properties Lagrangian Relaxed Assignment Problem \ .] 

For a function evaluation of & 29 we show that an optimal (or 

suboptimal) solution of this relaxed problem ([5] 6.2) can be 
constructed from that of a two-dimensional assignment problem. 
Then, the nonsmooth characteristics of @ 3 are addressed, 

15 followed by a method for computing a subgradient . 



Evaluation of 0 3 . Define for each {i 1 ,i 2 ) 

2 

j 3 = j 3 (i x , i 2 ) and a new cost function c ±i by: 



an index 



j 3 (i lf i 2 ) =arg min^c^^u^ \ i 3 = 0, 1 7 M 3 J, 
<i 2 = <i 2 i 3 (V 2 ) + < for Ux.ijJMOrO), 



( [5] 6.3) 
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Given an index pair (i x ,i 2 ) , j 3 need not be unique, resulting 
in the potential generation of several subgradients ([5] 6. 11) . 
(We further discuss this issue at the end of the Subsection [ 
5.2.3] . ) 
5 Now, define 

<P 3 (u 3 ) = Minimize 0 3 (z 2 ;u 3 ) = £ £ cfizf,-^ uf. 

i 1 =0 i 2 =0 i 3 =0 

Subject To ^z±±=l, i 1 = l f ... f M 1 , (6.4) 

I 2 =0 1 2 

4 s 2 ■ - 

z? , e { 0 , 1 } for all i, , i 9 . 

12 1 ^ 



[ (5.4) ] 

As an aside, two observations are in order. The first is that 
the search procedure needed for the computation of the relaxed 
cost coefficients in ([5] 6.3) is the most computationally 
10 intensive part of the entire relaxation algorithm. The second 
is that a feasible solution [z 3 ]zi of a sparse problem ([5]6..1) 
yields a feasible solution [z 2 ]z£ of ([5] 6.4) via the 
construction 



2 

z. 



1/ if zl^.i^ 1 for some (i lf i 2 ,i 3 ), 
0, otherwise. 
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Thus the two-dimensional assignment problem ( [5] 6.4) has 
feasible solutions other than those with exactly one nonzero 
index if the original problem ([5] 6.1) does. 

The following Theorem [5] 6.1 states that an optimal 
solution of ([5] 6.2) can be computed from that of ([5] 6.4). 
Furthermore, if the solution of either of these two problems is 
e -optimal, then so is the other. The converse is contained in 
Theorem [5] 6.2. 

Theorem [5] 6.1. Let [co 2 ]^ be a feasible solution to the two^ 
dimensional assignment problem ([5] 6.4) and define w 3 by 



([5] 6.5) 



if i 3 =j 3 and (i lf i 2 ) * (0, 0) / 



w. 




if i 3 *j 3 and i 2 ) * (0,0) , 




r 00i 3 



= 0 



if Cq^ + u^O, 
if c 0 3 0i3 + u, 3 3 >0. 



Then w 3 is a feasible solution of the Lagrangian relaxed problem 
(6.2) and [03 (w>,u 3 ) = <p3 (w 2 ;u 3 ) ] <p 3 (w 3 ; u 3 ) = <p 3 {w 2 ;u 3 ) . If, in 
addition, w 2 is optimal for the two [ ] ^dimensional problem, then 
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w 3 is an optimal solution of the relaxed problem and 
$ 3 (u 3 ) =<P 3 (u 3 ) . 

With_ the exception,, of one equality being converted to 

5 an inequality, the _following theorem _is a_ converse of 

Theorem [5] 6.1. 

Theorem [5] 6.2. Let [w 2 ]^ be a feasible solution to problem 
([5] 6.2) and define w 2 by 

M ([5]6. [7]6) 

<i 2 = £ <i 2 i 3 for (1^^)^(0,0), 

i 3 =0 

^ 0 2 0 = 0 if (i u i 2 ) = (0,0) and + u^O for all i 3/ 
w£ 0 = l if (i lf i 2 ) = (0, 0) and Cq 0 ^ + u^z0 for some i 3 . 



Then w 2 is a feasible solution _of the proJblem ( [5] 6.4) and 

0 3 ( w 3 ; u 3 ) * <p 3 (w 2 ; u 3 ) . If, in addition, w 3 is optimal for 
([5] 6.2), then w 2 is an optimal solution of ([5] 6.4), 
<p 3 (w 3 ;u 3 ) =<p 3 (w 2 ;u 3 ) and # 5 ( u 3 ) = <£, ( u 3 ) . 

15 

[The Nonsmooth Optimization Problem.] 
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An Algorithm for Computing 0 3 and a Subgradient. 
Given [u 3 ]uf the problem of computing 0 3 ([u 3 ]uf) and a 
subgradient of 0 3 ([u 3 ]ui) can be summarized as follows. 

1. [ ] Starting with problem ( [5] 6.1) , form the relaxed problem 

( [5]6.2) [ U 

2. To solve ([5] 6.2) optimally, define the two [ ] ^dimensional 
assignment problem ([5] 6.4) via the transformation 
( [5]6.3) [ U 

3. Solve the two-dimensional problem ([5] 6.4) optimally. 

4. Reconstruct the optimal solution, say [w] t3] # 3 , of ([5] 6.2) 
via equation (6 . 5 f . 61 ) as in Theorem [5] 6.1. 

5 . Then 



* 3 (u 3 ) S ^(^;u 3 ) = £ f; f; o\ xi%i e w ^ 



i^O i 2 =0 i 3 = 0 



% < £ f: ^w,-! 



(6.7) 



[(5.10)] 

A subgradient is given by substituting [ w*]_Q 3 into the 
objective function ( [5] 6.2) , erasing the minimum sign, and 
taking the partial with respect to u? . The result is 
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3 

<3i = 



d0 3 {u 3 ) 



ttK 

jL=0 i,=0 



-1 



for i^ = l f 



(6.8) 



[ (5.11)] 
VI . 3 . The Recovery Procedure [ . ] 
The next objective is to explain a recovery procedure, 
5 i.e., given a feasible (optimal or suboptimal) solution w 2 of 
([5] 6.4) (or w 3 of ([5] 6.2) constructed in Theorem [5] 6.1), 
generate a feasible solution [z 3 ]zf of equation ([5] .6.1) which 
is close to w 2 in a sense to be specified. We first assume that 
no variables in ([5].6. 1) are preassigned to zero; this 

10 assumption will be removed shortly. The difficulty with the 
solution w 3 in Theorem [5]j>.l is that it need not satisfy the 
third set of constraints in ([5]6..1). (Note, however, that if 
w 2 is an optimal solution for ([5]6.4) and [ ] w 3 , as constructed 
in Theorem [5] 6.1, satisfies the relaxed constraints, then w 3 

15 is optimal for ( [5] 6..1) .) The recovery procedure described here 
is designed to preserve the 0-1 character of the solution w 2 of 
([5]J5.4) as far as possible: If i = 1 and 1^0 or i 2 *0 , the 
corresponding feasible solution z 3 of ([5] 6.1) is constructed 
so that zli^l for some (i 1 ,i 2 ,i 3 ). By this reasoning, 

20 variables of the form z 0 3 0l3 [,] can be assigned a value of one in 
the recovery problem only if w$ 0 - 1 . However, variables Zooi 3 
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will be treated differently in the recovery procedure in that 

they can be assigned either_ zero or one independent of the 

value Wq 0 . This increases the feasible set of the recovery 

problem, leading to a potentially better solution. 

5 Let (I 2 ) , i 2 ( 2 2 ) j^ 1 ^ be an enumeration of indices {i l7 i 2 } 

of w 2 (or the first 2 indices of w* constructed in equation 

(6^5) ) such that Wi 2 l(l2 >, i2 (i 2 > =1 for (i x {l 2 ) , i 2 d 2 ))* (0,0) and 

(iiUg) /i 2 (I 2 ))= (0,0) for [I 2 =0]1 2 ^0 regardless of whether W q 0 = 1 

or not. (The latter can only improve the quality of the 

10 feasible solution. ) 

Next to define the two-dimensional assignment problem that 

restores feasibility, let 

2 3 [ ( [5]6. [12] 9)] 

c i 2 ± 3 = c i 1 ii 2 )i 2 U 2 )i i for 12 = 0,...,!^. 



Then the two-dimensional assignment problem that restores 
15 feasibility is 



Minimize Jj Jj C L *L 



i,=o i,»o 2 3 2 3 



\S 2 



i 3 =o ^ 3 



Su£>jec 

£ z i,i =lf i 3 = l,...,M,, 



i2=0 « • 3 3 ' 

z^efo,!} foralll 2 ,i 



2' -3' 



( [5]6.1[4] 0) 
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The next objective is to remove the assumption that all 
cost coefficients are defined and all zero-one variables are 
free to be assigned. We first note that the above construction 
of a reduced r dimension! two-dimensional assignment problem 
5 ([5]6.1[3]1) is valid as long as all cost coefficients c 2 are 
defined and all zero-one variables in z 2 are free to be 
assigned. Modifications are necessary for sparse problems . If 
the cost coefficient c 2 ? <i >i (j ji is well-defined and the zero- 
one variable z% { X )± {1 )± is free to be assigned to zero or one, 

10 then define = c^,^,^ as in < [5] 6. [12] 9) with z^^^^^ 

being free to be assigned. Otherwise, z \t is preassigned to 
zero or the corresponding arc is not allowed in the feasible 
set of arcs. To ensure that a feasible solution exists, we now 
only need ensure that the variables z^ 0 are free to be assigned 

15 for i 2 = 0, 1, 1^ with the corresponding cost coefficients being 
well-defined. [(] Recall that all variables of the form z* 00 , z^^ 
and z 0 3 0i3 are free (to be assigned) and all coefficients of the 

3 3 3 

corresponding form c i00 , c oi ^ 0 and c Q01 ^ need to be defined . [)] 
This is acconplished as follows: [ If] if the cost coefficient ( 2 }1 (I )0 
20 is well-defined and z i x (i 2 )i 2 {i 2 )o ^ s free, then define 
c i 2 o = °i 2 (i 2 )i 2 [i 2 )o with z^q being free. Otherwise, since all 
variables of the form z* 0 and z^^ are free to be assigned 
with the corresponding costs being well-defined, set 
c* Q = c? (I)00 + Cq^^^q^ where the first term is omitted if i 1 (I 2 )=0 



231 



NUMERI .01 USC 2 

and the second, if i 2 (I 2 )=0. For i 1 (l 2 )=0 and i 2 (l 2 )=0 , define 

2 _ 3 
c oo ~ c ooo * 

VI . 4 . The final recovery problem \ . 1 
The recovery problem for the case riV=3l N = 3 is 

r <[5]6.1[4]1) 



Minimize £ £ c^z^ 



i 2 -o i 3 =o ^ 3 ^ 3 



Subject to 5^ z^i =1/ 1 2 = 1,...,L 2 , 



2 z i 2 2 i 3 =1 ' . 

2 2 =0 2 3 

z*^ 6 {0,1} for all l 2 ,i 3 . 



1=1,...,AL, 



Let Y be an optimal or feasible solution to this [2- 
dimensionall two-dimensional assignment problem. The recovered 
feasible solution z 3 is defined by 

( [5]6.1[5]2) 



,3 J 1 ' if " WV' W 1 : 
Ji iVs \o, otherwise. 



) and Y. =1, 



10 
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Said in another way, let {(2 2 (I 3 ) , i 3 f^)}* 3 ^ be an enumeration 
of indices {l 2 ,i 3 } of Y such that Y i 2 a 3 ) ,i 3 u 3 ) ^ or 
(1 2 (1 3 ) , i 3 (l 3 )l*(0,0) and (I 2 (I 3 ) , i 3 (I 3 1) = (0, 0) for I 3 =0 regardless 
of whether Yqo- 1 or not. Then the recovered feasible solution 
5 can be written as 

«i 3 Ii={j' ^ V^^VV^VV-^ ( [5 ]6.1[6]3) 

[ 0 , o therwise . 



10 



The upper bound on the feasible solution is given by 

V 3 (2 3 )=f2 fl £ <V 3 <i 2 i3 (t5]6.1[7]4) 

i 1= 0 i 2 = 0 i 3 = o 123 

and the lower by & 3 ([u 2 ]ui) for any multiplier value ( [u 3 ] uf) . 
In particular, we have 
15 & 3 (u 3 ) <. & 3 {u 3 ) < V 3 (z* 3 ) <. V 3 (2 3 ) , where u 3 is any multiplier 

value, [(u 3 )]0^ is a maximizer of 0 3 ([u 3 ]uf), —3 is an optimal 
solution of the multidimensional assignment problem and ^ 2 3 
is the recovered feasible solution defined by ( [5] 6.. [16) . 

Other Relaxations for the Multidimensional Assignment 
20 Problem, (ABP) ] 13) . 

VII. OTHER RELAXATIONS FOR THE MULT I D I MENT I ONAL 

ASSIGNMENT PROBLEM 

In this section, we briefly describe other possible 
relaxations and their implications. These include linear 
25 programming relaxations and the corresponding duals. 
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Recall from Section II that one starts with the definition 
of the zero-one variable z. . =1 and cost coefficients to form 

X tT X lt 

the problem [ @ ] 

Minimize \\ c. . z. . 

■ L l" ± N 

Subject To 2 z l i =1 (i,=l,-»^)f 

J =1 (i 2 = l,~,M 2 ) , (7.1) 

ilia™ 1 !* 1 " 



10 



(i p =l,...,M p and p = 2,-,N-l), 



15 w h 

ere [c 0 . 0 ]c 0 0 is arbitrarily defined to be zero. Here, each 
group of sums in the constraints represents the fact that each 
non-dummy report occurs exactly once in a "track of data". One 
can modify this formulation to include multiassignments of one, 

20 some, or all the actual reports. The assignment problem 
([2J3..5) is changed accordingly. For example, if z i ^ is to be 
assigned no more than, exactly, or no less than n ijt times, then 
the " = 1" in the constraint ([2]3.5) is changed to "£,=,*n*, lf 
respectively. Modifications for group tracking and 

25 multiresolution features of the surveillance region will be 
addressed in future work. In making these changes, one must 
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pay careful attention to the independence assumptions that need 
not be valid in many applications. 

Again, the recent work of Poore and Drummond (A. B . Poore 
and O.E. Drummond. Track initiation and maintenance using 
5 multidimentional assignment problems. In D.W. Hearn, W.W. 
Hager, and P.M. Pardalos, editors , Network Optimization, volume 
450, pages 407-422, Boston, 1996. Kluwer Academic Publishers 
B . V. ) significantly extends the formulation of the track 
maintenance and initiation to new approaches. The discussions 
10 of this section apply equally to those formulations. 

n VII . 1 . The Linear Programming Relaxation \ . 1 

In the linear programming relaxation, one replaces the 

zero-one variables constraints 

z. . 6(0,1} for all i,... f i„ ([6]7.2) 
15 with the constraint 

Osz. . si for all i lf ...i N . ([6] 7. [2] 3) 

Then, the problem ([6] 7.1) can be formulated as a linear 
programming problem with the constraint ([6] 7.2) in ([6]7.-l) 
replaced by ([6]7.3) with a special block structure to which 
20 the Dantzing-Wolf e decomposition is applicable. Of course, 
after solving this problem, one must now recover the zero-one 
character of the variables in ( [6] 7.1) and there are many ways 
to do this, such as using the two [ ] ^dimensional assignment 
problems. Commercial software is also available. 

25 
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VII . 2 . fl The Linear Programming Dual and 

Partial Lagrangian Relaxations \ . 1 

Given the linear programming relaxation, one can formulate 
the dual problem or the partial Lagrangian relaxation duals 
5 with respect to any number of constraints. In particular, this 
is precisely what is done in Section 3 on the Lagrangian 
relaxation algorithm presented. The much broader class of 
algorithms provided in the US Patent 5537119 (Aubrey B. Poore, 
Jr., Method and System for Tracking Multiple Regional Objects 

10 by Multidimensional Relaxation, US Patent Number 5537119, 
filed 14 March 1995, issue date of 16 July 1996) can be 
modified to remove the zero-one character when one relaxes M 
sets of constraints to an ±N-M[ dimensional 1 ) -dimensional 
problem and recovers with an [N-M^-l dimensional! (N-M +1) - 

15 dimensional problem. This avoids the problems associated with 
the NP-hard (N-M) - and fN-M+11 (N-M +1) - dimensional problems. 
However, to restore the zero-one character, one can do it 
sequentially with an assignment problem or with one of the many 
zero-one rounding techniques. These formulations are easy to 

20 work out and thus the details are omitted. 

The complete dual problem is another way of solving the LP 
problem and may indeed be more efficient in certain 
applications. In addition, the solution of this dual problem 
may provide excellent initial approximations to the multipliers 

25 for Lagrangian relaxations. 
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[Multisensor Resolutions Problems and Their Solution via 
Linear Nonlinear Assignment Problems and Inequality 
Constraints. (ABP) The mathematical development for 
inequality constraints is understood by Aubrey B. Poore, but 
5 this work has not been written up. This feature is important 
in combining information from sensors with different 
resolutions such an IR and Radar. 



Formulation and Algorithms for the Distributed Sensor 
10 Tracking: Track- to-Track Correlations. (ABP) Hot Starts 

for Track Maintenance and Initiation: Bundle Modifications 
(ABP) .] 

VIII. HOT STARTS FOR TRACK MAINTENANCE 

AND INITIATION: BUNDLE MODIFICATIONS 

15 Thus reuse of multipliers and the first proof that this 

reuse is actually valid was presented in this section. The 
reuse in track maintenance is demonstrated and discussed in the 
work of Poore and Drummond [ [xx] ] (A.B. Poore and O.E. Drummond. 
Track initiation and maintenance using multidimensional 

20 assignment problems. In D.W. Hearn, W.W. Hager, and P.M. 
Pardalos, editors, Network Optimization, volume 450, pages 407- 
422, Boston, 1996. Kluwer Academic Publishers B.V.) . The only 
item left is the modification of the bundle of subgradients for 
the use with the multiplier values as one goes through the 

25 recovery problem as in Section IV. 3 . 6 or as one moves the 
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window in track maintenance. It is this aspect of the 
nonsmooth optimization that adds an order of magnitude to the 
speed of the algorithms. This work is based on both a 
mathematical proof as in Section IV. 3 . 6 as well as 
5 computational experience and heuristics . 

IX . [Adaptive Window Size Selection 

(ABP) « 1 ADAPTIVE AND VARIABLE WINDOW SIZE SELECTION 

These data association algorithms, based on the 
multidimensional assignment problem, range from sequential 

10 processing to many frames of data processed all at once! The 
data association problem for multisensor[ ] nmiltitarget 
tracking is formulated using a window sliding over the frames 
of data from the sensors. This discussion focuses on the work 
of Poore and Drummond and not on the formulations in Section 

15 [2] III. Firm data association decisions for existing track are 
made at, say frame k, with the most recent frame being Jc+iVT. 
Decisions after frame k are soft decisions. Reports not in the 
confirmed tracks are used to initiate tracks over frames 
numbered k-I to 7c+M. 

20 The length of these windows varies from sequential 

processing, which corresponds precisely to [1=0] I = 0 and 
\M=1] M = 1 , to user defined large values of J and [ . (In the 
case of sequential processing, we also have a temporary track 
file of dynamically feasible tracks, but incorrect data 

25 association.)] There is a marked change in performance over 
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this range. As the two windows become exceedingly long, there 
is little statistically significant information gained and 
indeed performance can degenerate slightly. Thus, one must 
manually find the correct window lengths for performance in a 
5 given scenario and the algorithms do not change for a changing 
environment. Thus we propose an adaptive method for adjusting 
the window lengths. (The method has been highly successful for 
selecting the correct window length for multistep methods in 
ordinary differential equations.) 
10 1. Compute the solution for different window lengths that 
overlap and differ by one or two frames. 



2. Compare the solution quality, i.e., the value of the 
objective function, for two adjacent window lengths. 

15 

3. If there is a predefined improvement in either direction, 
we then, for stability, repeat the exercise for a shorter 
or longer on the first try. If there is consistent 
information, we adjust the window size in the indicated 

20 direction. This can be given a well defined mathematical 

formula in terms of the assignment problems of different 
dimensions . 

[• Algorithm Enhancements due to Data Structures (ABP) . 
• Construction and enumeration of the Best K-Solutions. 
25 • New Nonsmooth Optimization Techniques. 
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SENSITIVITY ANALYSIS 

For sequential processing in which the two [ ] ^dimensional 
5 assignment algorithm is used, one can use the LP sensitivity 
analysis theory and easily obtain the corresponding answers. 
For the multiframe processing, the optimal solution is not 
available; however, there are still two approaches one can use. 
First, the basic algorithm can perform the same sensitivity 

10 analysis at each stage (loop) of the algorithm as is done in 
the two [ ] ^dimensional assignment problem since the evaluation 
of function ® N . k+1 is equivalent to a two [ ] ^dimensional 
assignment problem. Alternately, one could use an LP 
relaxation of the assignment problem and base the sensitivity 

15 on the resulting LP problem. We currently see this as an 
important step in finding even better solutions to the 
assignment problem if so desired. 

XI. [New Auction Algorithms (ABP & PS) 1 NEW AUCTION 

20 ALGORITHMS . 

In this section we present a new auction algorithm for the 
two [ ] ^dimensional assignment problem. 

An important step in solving N-dimensional assignment 
problems for N > 3 [ , ] is finding the optimal solution of the 
25 [2 -dimensional 1 two-dimensional assignment problem. In 
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particular, we wish to solve the r2-dimensionall two-dimensional 
problem which includes the zero -index. T [3 , ] ypically this 
problem can be thought of as trying to find either a minimum 
cost or maximum utility in assigning person to objects, tenants 
5 to apartments, or teachers to classrooms. We will follow the 
work of Bertsekas and call the first index set persons and the 
second index set objects. The statement of this problem is 
given below ( 11. 1) when the number of persons is m and the 
number of objects is n. 

10 I (11.1) 

Minimize £ £ c n K n 
i=o j=0 J 

Subject to ^2 x.. = l for all i = l,...,m, 



x. . = 1 for all j = 1 , . • . , n> 

i = 0 1J 



There are a couple of assumptions which we will make about 
( 11 . 1) . First of all, we assume that [c i0 ] c i0 and [c 0 ^] c 0j - are 
well-defined and the corresponding variables t^io^io anc * 
[x 0 j]x Q j are free to be assigned. Second if a cost c ± j happens 
15 to be undefined, then the corresponding variable [x^x^ will 
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be set to 0. In effect if [Cylc^ is undefined, we simply 
remove this cost and variable from inclusion in the problem. 

Notice that there are no constraints on the number of 
times person 0 and object 0 can be assigned. But notice that 
5 the first constraint requires that each non-zero person i must 
be assigned exactly once. Similarly, the second constraint 
forces each non-zero object j to be assigned exactly one time. 
Finally, the last constraint gives an integer solution, 
although we will see shortly that this constraint can be 
10 relaxed to admit any solution rx^^Ql x^Q . One reason for 
requiring that all of the costs of the form [c i0 ]c i0 and [c 0 j]c 0 j 
be defined is so that we are guaranteed a feasible solution 
exists for the given problem. 



XI . 1 . fl Relaxation of One Constraint f . 1 



15 



Relaxing the last constraint, we obtain the following 



problem: 



Minimize 



i=0 j=0 





Subject to 




for all i = 1, 



V{o # i). 
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This lower bound is achieved and the second constraint in 
the original problem is satisfied if the following conditions 
are met : 



For i=0, 



For i*0, 



X = i 



1, if c w + u^0 
0, otherwise 



1, if j = arg min{c ik + u£\k = 0,l, n } 
0 otherwise 



All of j f s are assigned only one time. 



5 The relaxation of the first constraint is analogous and would 
lead to similar results with i and j exchanged and the 
introduction of the multiplier u 1 . 

XI . 2 . [] Relaxation of Both Constraints T . 1 

This time we will relax both constraints and [using! use duality 
10 theor y to obtain a relationship between the multipliers u 1 and 
u 2 . 

[ a bunch of equations, theorems and the such here. 

From the above relaxation, we arrive at the following 
complementary slackness conditions which have been relaxed by 
15 a small value 6.] 
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Definition: An assignment S and a multiplier set (u 1 ^ 2 ) are 
said to satisfy e- complementary slackness (e - CS) if 



c if + 1^ + 1^ = 0 for all (i,j)eS, 
c ±j + ul + Uj>-e for all (i,j)eA. 

Forward Auction Algorithm [ . ] 

(1) Select any unassigned person i 
5 (2) Determine the following quantities: 

j. = arg min j c. k + u£\k eA(i) J, 

2 

w.= min \ c ik + ul\k 6A(i) , 



In the selection of j ± above, if a tie occurs between 0 
and any non-zero index k, then select j L as 7c. Otherwise, 
if there is a tie between two or more non-zero indices, 
10 the choice of j ± is arbitrary. 

Also if A(i) consists of only one element, then set w i =°°. 

(3) Update the multipliers and the assignment: 
If Ji = 0, then 

(a) Add (1,0) to the assignment. 
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(b) Update u*: = -c i0 . 
If jfi * 0, then 

(a) Add (i,j ± ) to the assignment. 

(b) Remove (i', j^) from the assignment if j\ was 
5 previously assigned. 

( c ) Update u £ : = uj\ + ( w ± - v ± ) + e = w ± - c ±j + e . 

(d) Update := c^. + = -ft^-e. 



10 Reverse Auction Algorithm. 

'0j 



(1) Select any unassigned object j, such that c Q . + u ? >0< 



(2) Determine the following quantities: 
i. = arg min j c jk + u£\k eB{j) J, 



In the selection of i ^ [ . ] above, if a tie occurs between 
15 0 and any non-zero index k, then select j i as k. 

Otherwise, if there is a tie between two or more non-zero 
indices, the choice of j ± is arbitrary. 
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Also if B(j) consists of only one element, then set 
Yj = 00 . 

(3) Update the multipliers and the assignment: 
If ij = 0, then 

(a) Add (0,j) to the assignment. 

(b) Update u *: = -c ojm 
If ij * 0, then 

(a) Add to the assignment. 

(b) Remove from the assignment if was 
previously assigned. 

(c) Update u*. : = u*. + (y. - (3..) + e = y. - c ±j + e . 

(d) Update u 2 : = - ( c. ^ + u^) = -y^e . 
Combined Forward/Reverse Auction Algorithm 

1. Assume that u 2 is given as an arbitrary multiplier. 

2. Adjust the value of u 2 for each object j as follows: 

If c 0j + <0, then set u 2 : = -c Qj . - 

3. Run iterations of the Forward Auction Algorithm until 
all persons become assigned. 

4. Run iterations of the Reverse Auction Algorithm until 
all of the objects become assigned. 

Note at the completion of the Forward auction step we have 
the following conditions satisfied: 

• c 0j + Uj ^0 for all objects j. 

• c-. + Ui+u 2 =0 for all S. 
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• c ±j + u j * min { c i* + u k | + e for all (i,j)e S. 
Thus we can prove the following proposition. 
Proposition : If we assume that c 0 _. + u? ^ 0 at the start of the 
Forward Auction Algorithm and all of the persons are assigned 
5 via a forward step, then we have: 

c^. + u^ + u* k -6 for all A. 

c ± . + ul + Uj =0 for all (i,j)e 5. 

c^ + uf <; min |c iic + J + e for all S. 

10 Optimality of the Algorithm 

Theorem: e-CS preserved during every forward and reverse 
iteration. 

Theorem: If a feasible solution exists, then the resulting 
solution is with me of being optimal for the Combined 
15 Forward Reverse Algorithm. 

[ Implementation Specifics 
Parallelization. 

Here are but a few comments.] 
20 XII. SOME CONCLUDING COMMENTS 

Although the algorithm appears to be serial in nature, its 
primary computational requirements are almost entirely 
parallelizable . Thus parallelization is planned. 
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Step 2 is the computational part of the algorithm. 
Evaluating ® N _ k+1 and computing a subgradient use in the search 
procedure requires 99% of the computing time in the algorithm. 
This part uses two [ ] ^dimensional assignment algorithms, a 
5 search over a large number of indices, and a nonsmooth 
optimization algorithm. It is the second part (the search) 
that consumes 99% of the computational time and this is almost 
entirely parallelizable . Indeed, there are two [ ] ^dimensional 
assignment solvers that are highly parallelizable. Thus, we 

10 need but parallelize the nonsmooth optimization solver to, have 
a reasonably complete parallelization. 

If a sensitivity analysis is desired or if one is 
interested in computing several near-optimal solutions, a 
parallel processor with a few powerful processors and good 

15 communication such as on the Intel Paragon would be most 
beneficial . 

The foregoing discussion of the invention has been 
presented for purposes of illustration and description. 
Further, the description is not intended to limit the invention 

20 to the form disclosed herein. Consequently, variation and 
modification commensurate with the above teachings, within the 
skill and knowledge of the relevant art, are within the scope 
of the present invention. The embodiment described hereinabove 
is further intended to explain the best mode presently known of 

25 practicing the invention and to enable others skilled in the 
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art to utilize the invention as such, or in other embodiments, 
and with the various modifications required by their particular 
application or use of the invention. It is intended that the 
appended claims be construed to include alternative embodiments 
5 of the invention to the extent permitted by the prior art. 
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What is claimed is : 

1. A method for tracking a plurality of objects, 
comprising: 

repeatedly scanning a region containing a set consisting 
5 of one or more moving objects and generating N sequential 
images or data sets of said region, a plurality of observations 
in said images or data sets providing positional information 
for objects in said set; 

determining a plurality of tracks, at least one track for 
10 each object in said set; 

determining a plurality of costs, wherein each cost is for 
assigning one of said observations to one of said tracks; 

defining a linear programming problem: 

[ Minimize X^...^ • • i w z i l - * 

15 Subject To ^...i/i^'-i/ 1 (i x = l,...,*^) 

^V 3 ...i N Z ir • *i w = 1 (i 2 = 1, • • wM 2 ) 



/ ±± ± % i • • • i 1 

. .i p _ 1 i p+1 . . .i N i x i N 

(i p = l,...,M p and p=2, . . .N - 1) 

20 0 <> z i!...i N * 1 for a11 i lf ...,i N ,] 
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Minimize \\ c. . z. 

Subject To Y\ z, =1 (i. = 1, . . . , 



i i i i ' " " w 



(i 1 = l / . . . r M p and p = 2, . . . , N-l) , 



O^z. . ^ 1 for all 



wherein each c ± ± is included in said plurality of costs, each 
M if i = l, . . . ,N, being one of: (a) a number of observations in an 
i th image or data set of said N sequential images or data sets; 
(b) a sum of a number of tracks in said plurality of tracks, 
5 and a number of said observations in the i th image or data set 
not assigned to one of said tracks; and (c) a number of tracks 
in said plurality of tracks; 

solving said linear programming problem for values of 
z ± < for each i 2 ...i N ; 
10 determining a value z ± ±N in {0,1} for each 

corresponding to each z ± ± , wherein said values z ± t provide 
an optimal or near optimal solution to said linear programming 
problem; 
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taking one or more of the following actions based on said 
optimal or near-optimal assignment of said plurality of points 
to said plurality of tracks: 

[ ] sending a warning to aircraft or a ground or 

5 sea facility, 

controlling air traffic, 

controlling anti-aircraft or anti-missile equipment, 
taking evasive action, 

working on one of said one or more objects, 
10 surveilling one of said one or more objects. 
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